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Abstract 


The  micro- Doppler  (m-D)  effect  appears  in  the  synthetic  aperture  radar  (SAR)/inverse 
SAR  (ISAR)  image  of  a  target  whenever  the  target  has  one  or  more  rotating  or  vi¬ 
brating  parts.  M-D  effect  introduces  distortion  in  the  SAR/ISAR  images.  On  the 
other  hand,  m-D  effect  also  carries  information  about  the  features  of  moving  parts  of  a 
stationary  or  moving  target  that  can  be  used  for  target  identification  purpose.  Based 
on  L-statistics,  spectrogram,  and  inverse  Radon  transform,  a  radar  data  exploitation 
Matlab  toolbox  was  develped  for  classifying  air,  land  and  ocean  targets,  which  helps 
Geospatial  Intelligence  (GEOINT)  support  to  Canadian  Forces.  This  toolbox  will 
scan  and  search  any  M-D  activity  in  the  desired  area  of  the  SAR/ISAR  image.  The 
toolbox  has  three  options;  1)  Focus  the  stationary  target  if  there  is  no  m-D,  2)  Focus 
the  target  after  removing  the  m-D,  and  3)  extract  the  m-D  and  determine  the  motion 
parameters.  The  toolbox  was  tested  and  validated  against  a  variety  of  simulated  and 
measured  targets. 

Significance  for  defence  and  security 


The  micro-Doppler  (m-D)  effect  appears  in  the  synthetic  aperture  radar  (SAR)  image 
of  a  target  whenever  the  target  has  one  or  more  rotating  or  vibrating  parts.  Similar 
effect  appears  in  the  inverse  synthetic  aperture  radar  (ISAR)  imaging,  as  well.  If  the 
frequency  modulations  on  the  returned  signal  caused  by  the  moving  parts  are  not 
filtered,  then  the  m-D  effect  can  introduce  distortion  in  the  SAR/ISAR  images.  The 
frequency  content  of  the  m-D  signal  changes  over  time  in  a  wide  range.  Therefore, 
the  m-D  may  cover  the  rigid  body  and  make  it  difficult  to  detect.  On  the  other  hand, 
the  m-D  effect  also  carries  useful  information  about  the  features  of  moving  parts  of  a 
target  (type,  velocity,  size,  etc.)  that  are  complementary  to  existing  target  recognition 
methods.  It  is  easier  to  estimate  these  features  if  the  m-D  effect  is  separated  from 
the  rigid  body  part  of  the  radar  image.  Thus,  the  extraction  of  m-D  effects  from 
the  radar  images  have  attracted  significant  research  attention,  especially  the  image 
enhancement  and  m-D  extraction  for  the  spacebourse  sensors  such  as  RADARSAT-2 
and  TerraSAR-X. 

Based  on  L-statistics,  spectrogram,  and  inverse  Radon  transform,  a  radar  data  ex¬ 
ploitation  Matlab  toolbox  was  develped  for  classifying  air  and  land  targets.  This 
toolbox  will  scan  and  search  any  m-D  activity  in  the  desired  area  of  the  SAR/ISAR 
image.  The  toolbox  has  three  options;  1)  Focus  the  stationary  target  if  there  is  no  m- 
D,  2)  Focus  the  target  after  removing  the  m-D,  and  3)  extract  the  m-D  and  determine 
the  motion  parameters  such  as  rotation/vibration  rate,  initial  phase,  and  amplitude 
of  the  stationary  and  moving  targets.  The  toolbox  was  tested  and  validated  against 
a  variety  of  simulated  and  measured  targets. 


DRDC-RDDC-201 4-R6 


i 


The  m-D  toolbox  can  be  used  to  identify  important  classes  of  battlefield  targets. 
The  m-D  signatures  of  vibrating  and  rotating  structures  can  be  used  to  differenti¬ 
ate  and  identify  specific  types  of  military  vehicles/tanks/trucks  and  determine  their 
movement  and  the  speed  of  their  engines.  This  approach  could  provide  improved 
target  signature/feature  recognition  over  current  audio  classification  techniques.  In 
addition  to  battlefield  land  targets,  important  classes  of  air  and  ocean  targets  can 
also  be  distinguished  by  their  m-D  signatures,  for  example,  radar  signals  returned 
from  military  targets  that  incorporate  vibrating  or  rotating  structures  such  as  air¬ 
craft  propellers,  helicopter  rotors,  rotating  antennas  on  a  ship  or  an  aircraft  contains 
m-D  signatures  related  to  these  structures.  One  of  the  important  applications  of 
the  m-D  is  the  detection  of  a  hovering  helicopter  by  hub  radar  returns  through  m-D 
features  rather  than  the  blade  flash  returns,  which  is  currently  used  by  conventional 
scanning  radars.  Radar  returns  from  the  engine  compressor  and  blade  assemblies  of 
a  jet  aircraft  contain  jet  engine  modulation  that  also  generates  m-D.  Hence  the  m-D 
toolbox  offers  a  new  way  for  the  analysis  of  military  target  signatures.  It  provides 
additional  and  unique  target  features  that  are  complementary  to  those  made  available 
by  existing  methods.  The  toolbox  will  be  available  to  the  CFJIC  (Canadian  Forces 
Joint  Imagery  Centre)  and  MCE  (Mapping  and  Charting  Establishment)  clients  for 
Geospatial  Intelligence  (GEOINT)  support. 

Radar  gait  image  analysis  integrated  with  m-D  can  provide  new  identification  meth¬ 
ods  for  remote  detection  of  walking  personnel  either  in  battlefield  or  urban  scenarios. 
This  technique  can  be  applied  to  other  areas  as  well:  countering  terrorism,  conduct¬ 
ing  urban  military  operations,  providing  urban  border  security,  dealing  with  hostage 
acts,  intercepting  suicide  bombers,  and  detecting  humans  (soldiers)  in  a  forest.  Dis¬ 
tance  is  the  key  factor,  because  terrorists  and  other  criminals  often  act  from  afar, 
and  the  earlier  detection  means  a  quicker  response.  Therefore,  as  a  new  identi¬ 
fication/recognition  tool  for  operational  radar  Automatic  Gait  Recognition  (AGR) 
system,  m-D  has  the  good  potential  for  significantly  improving  military  capabilities 
and  protection  for  CF  and  CF  facilities. 


ii 
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Resume 


L’effet  micro-Doppler  (m-D)  apparait  dans  l’image  d’un  radar  a  synthese  d’ouverture 
(SAR)  ou  d’un  radar  a  synthese  d’ouverture  inverse  (ISAR)  d’une  cible  qui  comporte  au 
moins  une  piece  tournante  ou  vibrante.  L’effet  m-D  introduit  de  la  distorsion  dans  les 
images  SAR  ou  ISAR.  Par  contre,  cet  effet  foumit  aussi  des  renseignements  sur  les 
caracteristiques  des  pieces  mobiles  d’une  cible  fixe  ou  en  mouvement,  renseignements  qui 
peuvent  servir  a  identifier  la  cible.  A  partir  des  statistiques  lineaires,  du  spectrogramme  et 
de  la  transformee  inverse  de  Radon,  une  boite  a  outils  Matlab  d’exploitation  des  donnees 
radar  a  ete  mise  au  point  pour  classer  les  cibles  aeriennes,  terrestres  et  oceaniques  afin 
d’aider  au  soutien  du  renseignement  geospatial  (GEOINT)  des  Forces  canadiennes.  Cette 
boite  a  outils  balaie  la  zone  desiree  de  l’image  SAR  ou  ISAR  a  la  recherche  d’activites  m- 
D.  Elle  comporte  trois  options  :  1)  la  focabsation  de  la  cible  fixe  en  l’absence  de  m-D;  2) 
la  focabsation  de  la  cible  apres  1’ elimination  des  m-D;  3)  1’ extraction  des  m-D  et  la 
determination  des  parametres  de  mouvement.  Elle  a  en  outre  ete  mise  a  l’essai  et  validee 
en  fonction  de  diverses  cibles  simulees  et  mesurees. 


Importance  pour  la  defense  et  la  securite 


L’effet  micro-Doppler  (m-D)  apparait  dans  l’image  d’un  radar  a  synthese  d’ouverture 
(SAR)  d’une  cible  qui  comporte  au  moins  une  piece  tournante  ou  vibrante.  Un  effet 
similaire  apparait  aussi  dans  l’image  d’un  radar  a  synthese  d’ouverture  inverse  (ISAR).  Si 
les  modulations  de  frequence  sur  les  signaux  reflechis  causes  par  les  pieces  mobiles  ne 
sont  pas  filtrees,  l’effet  m-D  peut  introduire  des  distorsions  dans  les  images  SAR  ou 
ISAR.  De  plus,  le  contenu  frequentiel  du  signal  m-D  varie  en  fonction  du  temps  sur  une 
large  gamine  de  frequences.  En  consequence,  l’effet  m-D  peut  couvrir  un  corps  rigide  et  le 
rendre  difficile  a  detecter.  Par  contre,  il  contient  aussi  des  renseignements  utiles  sur  les 
caracteristiques  des  pieces  mobiles  d’une  cible  (type,  vitesse,  dimensions,  etc.), 
renseignements  qui  sont  complementaires  aux  methodes  existantes  de  reconnaissance  de 
cibles.  Ces  caracteristiques  sont  plus  faciles  a  evaluer  si  l’effet  m-D  est  separe  de  la  partie 
rigide  de  la  cible  dans  l’image  radar.  L’ extraction  des  effets  m-D  des  images  radar  fait 
done  l’objet  de  beaucoup  de  recherches,  surtout  relativement  a  1’ amelioration  des  images  et 
a  1’ extraction  des  m-D  pour  les  capteurs  spatiaux  comme  RADARS AT-2  et  TerraS AR-X. 

A  partir  des  statistiques  lineaires,  du  spectrogramme  et  de  la  transformee  inverse  de 
Radon,  une  boite  a  outils  Matlab  d’exploitation  des  donnees  radar  a  ete  mise  au  point  pour 
classer  les  cibles  aeriennes  et  terrestres.  Cette  boite  a  outils  balaie  la  zone  desiree  de 
l’image  SAR  ou  ISAR  a  la  recherche  d’activites  m-D.  Elle  comporte  trois  options  :  1)  la 
focabsation  de  la  cible  fixe  en  l’absence  de  m-D;  2)  la  focabsation  de  la  cible  apres  le 
retrait  des  m-D;  3)  l’extraction  des  m-D  et  la  determination  des  parametres  de  mouvement 
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comme  la  vitesse  de  rotation  et  de  vibration,  la  phase  initiale  et  l’amplitude  des  cibles 
fixes  et  mobiles.  Elle  a  en  outre  ete  mise  a  l’essai  et  validee  en  fonction  de  diverses  cibles 
simulees  et  mesurees.  La  boite  a  outils  m-D  peut  etre  utilisee  pour  identifier  des  classes 
importantes  de  cibles  sur  le  champ  de  bataille.  Les  signatures  m-D  de  structures  vibrantes 
et  toumantes  peuvent  servir  a  differencier  et  a  identifier  des  types  particuliers  de  vehicules, 
de  chars  et  de  camions  militaires, 

ainsi  qu’a  determiner  leur  mouvement  et  la  vitesse  de  leurs  moteurs.  Par  rapport  aux 
techniques  actuelles  de  classification  audio,  cette  methode  pourrait  ameliorer  la 
reconnaissance  de  la  signature  et  des  caracteristiques  d’une  cible.  En  plus  des  cibles 
terrestres  sur  les  champs  de  bataille,  des  classes  importantes  de  cibles  aeriennes  et 
oceaniques  peuvent  aussi  etre  distinguees  au  moyen  de  leurs  signatures  m-D.  Les  signaux 
radar  renvoyes  par  des  cibles  militaires  dotees  de  structures  vibrantes  et  toumantes, 
comme  les  helices  d’aeronef,  les  rotors  d’helicoptere  et  les  antennes  rotatives  d’un  navire 
ou  d’un  aeronef,  comportent  des  signatures  m-D  associees  a  ces  structures.  Une 
application  importante  de  l’effet  m-D  est  la  detection  d’un  helicoptere  en  vol  stationnaire 
par  les  echos  radar  du  moyeu  de  rotor  au  moyen  de  leurs  caracteristiques  m-D  plutot  que 
par  les  eclairs  de  pales  qu’utilisent  aujourd’hui  les  radars  a  balayage  traditionnels.  Les 
echos  radar  provenant  des  aubes  et  du  compresseur  de  moteur  d’un  avion  a  reaction 
conti ennent  une  modulation  qui  produit  aussi  des  effets  m-D.  Par  consequent,  la  boite  a 
outils  pour  m-D  offre  une  nouvelle  methode  d’ analyse  des  signatures  de  cibles  militaires. 
Elle  fournit  en  outre  des  caracteristiques  de  cibles  supplementaires  et  uniques  qui  sont 
complementaires  a  celles  produites  par  les  methodes  actuelles.  Elle  sera  offerte  aux  clients 
du  Centre  d’imagerie  interarmees  des  Forces  canadiennes  (CIIFC)  et  du  Service  de 
cartographic  (S  Carto)  aux  fins  de  soutien  du  renseignement  geospatial  (GEOINT). 


L’ integration  de  1’ analyse  de  la  demarche  dans  les  images  radar  et  de  l’effet  m-D  fournit  de 
nouvelles  methodes  d’identification  pour  la  detection  a  distance  de  marcheurs  dans  des 
scenarios  de  bataille  ou  des  scenarios  urbains.  Cette  technique  peut  aussi  s’appliquer  a 
d’autres  domaines  :  lutte  contre  le  terrorisme,  conduite  d’operations  militaires  en  zone 
urbaine,  prestation  de  services  de  securite  a  la  frontiere  en  zone  urbaine,  interventions 
dans  les  situations  de  prise  d’otages,  interception  de  bombes  humaines  et  detection  d’etres 
humains  (soldats)  en  foret.  La  distance  constitue  un  facteur  cle,  car  les  terroristes  et  les 
autres  criminels  agissent  souvent  de  loin,  et  une  detection  plus  precoce  permet  une  reponse 
plus  rapide.  Par  consequent,  en  tant  que  nouvel  outil  d’identification  et  de  reconnaissance 
pour  le  systeme  operationnel  de  reconnaissance  automatique  de  la  demarche  (AGR),  l’effet 
m-D  peut  ameliorer  de  fa9on  considerable  les  capacites  militaires  et  la  protection  pour  les 
Forces  canadiennes  et  leurs  installations. 


IV 
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1  Introduction 


If  a  target  or  any  structure  on  a  target  has  mechanical  vibration  or  rotation  in  addi¬ 
tion  to  its  bulk  translation,  it  might  induce  a  frequency  modulation  on  the  returned 
signal  that  generates  sidebands  about  the  target’s  Doppler  frequency  shift.  This  is 
called  the  micro- Doppler  (m-D)  effect  [1-10].  Radar  signals  returned  from  a  target 
that  incorporates  vibrating  or  rotating  structures,  such  as  propellers  of  a  fixed-wing 
aircraft,  rotors  of  a  helicopter,  or  the  engine  compressor  and  blade  assemblies  of  a  jet 
aircraft,  contain  m-D  characteristics  related  to  these  structures  [1,2, 5, 9].  The  m-D 
effect  enables  us  to  determine  the  dynamic  properties  of  the  target  and  it  offers  a  new 
approach  for  the  analysis  of  target  signatures. 

Micro-Doppler  effect  appears  in  the  SAR/ISAR  image  of  a  target  whenever  the  target 
has  one  or  more  rotating  or  vibrating  parts.  If  the  frequency  modulations  on  the 
returned  signal  caused  by  the  moving  parts  are  not  filtered,  then  the  m-D  effect  can 
introduce  distortion  in  the  SAR/ISAR  images  [1,2, 5, 9].  The  observation  of  very  large 
distortions  from  experimental  SAR/ISAR  data  has  been  reported.  On  the  other  hand, 
the  m-D  effect  also  carries  information  about  the  features  of  moving  parts  of  a  target 
that  are  complementary  to  existing  target  recognition  methods.  Several  papers  have 
been  written  about  the  ways  to  deal  with  the  m-D  effect.  The  wavelet  analysis  of  the 
helicopter  and  human  data,  along  with  the  time-frequency  (TF)  representation  based 
imaging  system,  is  presented  in  [1,2].  Details  on  the  physics  of  the  m-D  effect,  with 
some  typical  examples,  are  given  in  [1-9].  A  method  for  separating  the  m-D  effect 
from  the  radar  image,  based  on  the  chirplet  transform,  is  proposed  in  [2, 10].  It  is 
easier  to  estimate  these  features  if  the  m-D  effect  is  separated  from  the  rigid  body 
part  of  the  radar  image.  Thus,  the  extraction  of  m-D  effects  is  an  important  problem 
in  the  radar  imaging.  Since  the  spectral  content  of  the  signal  corresponding  to  m-D 
effects  is  time- varying,  time- frequency  (TF)  analysis  techniques  are  an  appropriate 
tool  for  the  m-D  features  analysis  and  extraction  [1-28]. 

The  most  common  and  simplest  form  of  the  signal  that  corresponds  to  the  m-D  effect 
is  a  sinusoidal  frequency-modulated  (FM)  signal.  This  kind  of  the  m-D  signatures  can 
be  extracted  based  on  the  inverse  Radon  transform  of  its  TFR  [5].  This  transform 
projects  a  two-dimensional  (2D)  sinusoidal  signature  into  a  single  point  whose  position 
corresponds  to  the  parameters  of  the  sinusoidal  signature.  The  technique  based  on  the 
inverse  Radon  transform  can  be  used  for  other  nonsinusoidal  periodic  m-D  forms.  In 
general,  the  m-D  effect,  as  a  nonstationary  signal  with  time-varying  spectral  content, 
can  be  separated  from  the  rigid  body  by  using  the  L-statistics  on  the  TF  signal 
representations  [5].  The  L-statistics  based  separation  of  the  m-D  effects  from  the 
rigid  body  is  very  efficient  and  simple  for  any  m-D  form.  The  Viterbi  algorithm  can 
improve  the  m-D  tracking  in  high  noise  environments  [4], 
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In  Section  2,  a  method  for  the  estimation  of  m-D  motion  parameters  is  provided. 
In  this  section,  the  inverse  Radon  transform  is  reviewed  and  the  methodology  is 
provided  for  m-D  motion  parameters  estimation.  Its  application  to  the  cases,  with  a 
fraction  of  period  of  the  m-D  signal  being  available,  is  presented.  Since,  the  m-D  is 
not  necessary  of  the  form  of  a  sinusoidal  frequency  modulated  (FM)  signal,  a  simple 
algorithm  for  period  estimation  of  a  general  form  of  the  periodic  signals  is  presented 
in  this  Section,  as  well.  The  m-D  must  be  separated  from  a  rigid  body  in  order  to 
be  further  analyzed,  thus  the  L-statistics  based  algorithm  for  the  m-D  separation  is 
presented  in  Section  3.  The  m-D  toolbox  based  on  L-statistics,  spectrogram,  and 
inverse  Radon  transform  is  given  in  Section  4.  This  toolbox  will  scan  and  search  any 
m-D  activity  in  the  desired  area  of  the  SAR/ISAR  image.  The  toolbox  will  have  three 
options;  1)  Focus  the  stationary  target  if  there  is  no  m-D,  2)  Focus  the  target  after 
removing  the  m-D,  and  3)  extract  the  m-D  and  determine  the  motion  parameters 
such  as  rotation/vibration  rate,  initial  phase,  and  amplitude  of  the  stationary  and 
moving  targets.  The  toolbox  was  tested  and  validated  against  a  variety  of  simulated 
and  measured  targets.  In  Section  4,  the  time-frequency  toolbox  is  also  illustrated. 
Conclusions  are  given  in  Section  5. 

2  Estimation  of  micro-Doppler  motion 
parameters 


The  m-D  can  be  used  to  determine  properties  of  a  target  [1-28],  since  the  m-D  para¬ 
meters  are  related  to  the  parameters  of  the  corresponding  micro-motions.  Therefore, 
the  radar  m-D  signatures  are  of  a  great  potential  for  identifying  properties  of  un¬ 
known  targets.  The  most  common  form  of  the  m-D  is  a  form  of  sinusoidal  frequency 
modulated  (FM)  signal. 

Sinusoidally  modulated  signals  appear  in  many  applications,  like  in  radars  and  com¬ 
munications.  In  the  radar  signal  processing,  the  fast  rotating,  vibrating  or  oscillating 
parts  reflect  a  signal  causing  m-D  effect  in  a  form  of  sinusoidally  modulated  signal. 
In  practice  it  is  very  important  to  extract,  decompose,  and  estimate  parameters  of 
these  kinds  of  signals,  since  they  are  easily  related  to  the  physical  dimensions  and 
other  properties  of  the  moving  objects. 

In  this  section,  we  will  present  a  method  for  analysis  of  sinusoidally  modulated  com¬ 
ponents  based  on  the  inverse  Radon  transform  of  signal’s  time-frequency  represent¬ 
ation.  The  Radon  transform,  widely  used  in  computer  imaging  applications,  is  also 
used  in  time-frequency  for  projecting  Wigner  distribution  in  order  to  detect  linear 
frequency  modulated  signals  [33-37].  Here,  we  will  use  the  inverse  Radon  transform 
rather  than  the  Radon  transform.  Note,  that  the  behavior  of  the  direct  and  inverse 
Radon  transform  is  completely  different,  in  contrast,  for  example,  to  the  Fourier 
transform  [38].  Since  the  Radon  transform  of  a  two-dimensional  signal  containing 
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a  two-dimensional  delta  function  is  a  sinusoidal  pattern  with  amplitude  correspond¬ 
ing  to  the  distance  of  the  point  from  the  origin  and  the  initial  phase  corresponding 
to  the  phase  of  the  point  position,  then  it  is  obvious  that  a  sinusoidal  pattern  in 
the  time-frequency  plane  (produced  by  a  time-frequency  representation  of  sinusoid¬ 
ally  modulated  signal)  will  project  to  a  two-dimensional  delta  in  the  inverse  Radon 
transform.  This  is  obviously  an  optimal  transform  for  a  two-dimensional  sinusoidal 
pattern,  since  all  signal  energy  from  the  time-frequency  domain  will  be  projected  in  a 
single  point  in  the  inverse  Radon  transform  domain.  The  method  will  be  introduced 
on  monocomponent  sinusoidally  modulated  signal.  Then  it  will  be  extended  to  noisy 
and  multicomponent  signals  that  include  one  or  more  sinusoidal  patterns.  Finally  the 
method  will  be  applied  to  periodic  and  non-periodic  patterns  that  are  not  produced 
by  a  sinusoidally  modulated  signals  at  all.  The  examples  illustrate  the  efficiency  of 
the  presented  method. 

2.1  Radon  and  inverse  Radon  transform 

A  projection  of  a  two-dimensional  function  f(x,y )  onto  the  x-axis  is 

OO 

Rs(x)  =  j  f(x,y)dy.  (1) 

—  CO 

A  rotated  version  of  a  two-dimensional  signal  may  be  described  in  a  rotated  coordinate 
system,  by  a  coordinate  rotation  transform.  For  an  angle  a,  it  reads 


£  ' 

cos(a)  sin(o;) 

X 

C  . 

—  sin  (a)  cos(a) 

.  y . 

The  projection  of  a  function  f(x,y )  onto  £,  with  a  varying  rotation  angle  a,  is  the 
Radon  transform  of  the  signal  f(x,  y) 


f  f  «,CR 


f  f  /(A,  C)5(A  -  €)rfAdC- 


(2) 


Let  us  consider  simple  setup  where  the  analyzed  image  is  two-dimensional  delta  func¬ 
tion  located  at  the  point  f(x,y )  =  S(x  —  x0)5(y  —  y0)  in  x,y  domain.  Projection  of 
f(x,y )  onto  x  axis  is 


OO 

Rf(x,0)  =  J  f  (x,  y)dy  =  S(x 

— OO 


Xq). 
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For  an  arbitrary  direction  £0  =  xq  cos(a)  +  yo  sin(a),  C0  =  —xo  sin(o;)  +  yo  cos(a),  the 
function  /(£,  Q  =  6(£  —  Co)^(C  —  Co)  results  in  the  Radon  transform 

CO 

Rf((,a)=  f  =  W-(o) 

—co 

=  (5(C  -  (x0  cos  (a)  +  y0  sin(a))).  (3) 

Note  that  this  is  a  sinusoidal  pattern  in  a  two-dimensional  (£,  a)  domain,  with  the 
amplitude  \J +  y\  and  the  phase  ip  =  arctan(y0/.:c0).  Of  course,  the  Radon  trans¬ 
form  is  periodic  in  a  with  2tt.  Projections  for  0  <  a  <  n  are  sufficient  to  calculate 
all  transform  values.  By  knowing  all  the  projections,  for  0  <  a-  <  tt,  we  can  calculate 
the  two-dimensional  Fourier  transform  of  f(x,y).  It  means  that  we  can  reconstruct 
a  two-dimensional  function  f(x,y )  from  its  projections  or  integrals  (basic  theorem 
for  computed  tomography).  The  inverse  Radon  transform  may  be  calculated  in  the 
Fourier  domain  or  by  projecting  back  the  Radon  transform  (back- projection  method). 
Thus,  a  point  in  the  (x,  y)  domain  transforms  to  a  sinusoidal  pattern  in  Radon  trans¬ 
form  domain.  It  means  that  a  sinusoidal  pattern  will  be  transformed  into  a  point 
by  using  the  inverse  Radon  transform  (IRT).  When  all  energy  is  concentrated  into  a 
point,  then  its  parameters  estimation  is  very  robust  and  reliable. 

2.2  Parameters  estimation 

Let  us  now  consider  sinusoidally  frequency  modulated  signal 

A  ( ' . 

=  Ax  exp  I  j—  sm 

\  Jm 

It  is  known  that  a  time- frequency  representation  T(t,  c u)  of  a  given  signal  concentrate 
the  signal  energy  along  the  signal  instantaneous  frequency 

Wi{t)  =  2nAmcos(2nfmt  +  9m) 

i.e.,  this  signal  in  (t,  uS)  plane  is  presented  as  sinusoidal  pattern  image  (with  more  or 
less  deviations  depending  on  the  time-frequency  representation  used  to  transform  the 
signal  into  time- frequency  domain).  If  we  change  the  time  coordinate  with  tp  =  27i  frnt 
then,  from  the  previous  section  we  can  conclude  that  the  IRT  of  the  obtained  image 
T (<p/ (2n  fm) ,  co)  reduces  to  single  point  where  distance  from  origin  correspond  to 
modulation  parameter  Ami  and  the  angle  of  the  point  is  equal  to  9m-  In  this  way,  we 
can  accurately  estimate  modulation  parameters  Am  and  ()rn- 

The  modulation  parameter  fm  can  be  estimated  in  the  following  way.  Let  us  introduce 
coordinate  change  from  t  to  tp  as  tp  =  at  where  a  is  a  parameter.  Now  we  can  vary 


(27 Tfmt  +  9r 


(4) 
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the  parameter  a  within  some  range  of  possible  values  and  search  for  the  value  a  that 
produces  single  point  IRT.  In  that  case  we  know  that  a  =  2n fm  and  we  can  estimate 
the  modulation  parameter  fm.  In  this  procedure,  it  is  needed  to  find  the  value  of  a 
when  the  IRT  reduces  to  a  single  point  within  the  considered  domain.  This  mean 
that  the  IRT  is  ideally  concentrated  and  the  concentration  measures  [35]  can  be  used 
to  detect  that  we  reached  a.  The  range  for  a  should  be  wide  enough  to  include  2nfrn. 
Its  limits  could  be  determined  as  the  minimal  and  the  maximal  expected  2tt  fm  in  the 
considered  case. 

The  estimation  algorithm  is  summarized  as: 

Step  1.  Start  from  a  frequency  modulated  signal  x(t)  with  unknown  modulation 
parameters.  Assume  that  the  modulation  frequency  satisfies  /m;n  <  fm  <  /max, 
where  fmin  and  /max  are  constants. 

Step  2.  Calculate  the  time-frequency  representation  T(t,u)  of  x(t).  Here  we  can 
use  any  time-frequency  representation  [8]  concentrating  the  signal  energy  along  the 
instantaneous  frequency  in  the  time-frequency  plane.  The  result  of  this  step  is  a 
two-dimensional  time-frequency  image  of  the  considered  signal. 

Step  3.  Consider  a  set  of  possible  a  as  M  equally  spaced  values  between  27r/min  and 

271"  / 'max- 

Step  4.  For  each  a  within  the  considered  set,  introduce  coordinate  change  if  =  at 
and  calculate  the  IRT  of  the  image  T(f/a,  uj). 

Step  5.  Calculate  the  concentration  measure  /j  of  the  obtained  IRT  for  each  a  and 
find  a  that  provide  the  highest  concentration. 

Step  6.  Estimate  the  modulation  frequency  as  fm  =  a/(2n). 

Step  7.  Find  the  position  of  IRT  maximum  calculated  with  a ",  i.e.,  IRT  of  T(<p/ad, 
uS).  Denote  the  detected  coordinates  as  xmand  ym. 

Step  8.  Estimate  the  modulation  amplitude  as 

Am  =  x/Xm  +  Vm- 

Step  9.  Estimate  the  modulation  phase  as  9m  =  arctan  — . 

In  the  case  of  non-sinusoidally  modulated  signals,  producing  non-sinusoidal  patterns 
in  the  time- frequency  plane,  the  presented  approach  will  produce  the  closest  sinusoidal 
pattern  form,  as  it  will  be  shown  in  examples. 
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2.3  Method  implementation 

We  use  the  spectrogram  and  the  S-method  as  time-frequency  representations  in  the 
algorithm  Step  2.  The  spectrogram  is  defined  as  a  squared  modulus  of  the  short-time 
Fourier  transform.  In  the  discrete  time  domain,  it  reads 

SPEC(n,  k)  =  \STFT(n,  k)\2 

Nw 

STFT(n,  k)  =  Y ]  w{m)x(n  +  m)e~j^mk , 

m= 0 

where  w(n)  is  the  analysis  window  of  the  length  Nw.  Along  with  the  spectrogram, 
we  will  use  the  S-method  as  a  time  frequency  representation.  The  discrete  S-method 
is  of  the  form 


SM(n,  k )  =  \STFT(n,  k)\2 


T  2  Re 


'  L 

ST  FT  in,  k  +  p)STFT*(n,  k  -  p) 

.p= i 


) 


where  beside  the  time-domain  window,  used  in  the  STFT  calculation,  we  have  a  para¬ 
meter  L  that  corresponds  to  the  number  of  spectrogram  correcting  terms  [8].  It  is 
known  that  the  S-method  can  produce  highly  concentrated  time-frequency  represent¬ 
ation  of  a  given  signal.  The  S-method  is  numerically  very  efficient  since  there  is  no 
need  for  signal  oversampling.  Two  special  cases  of  the  S-method  are  the  spectrogram 
(with  L  =  0)  and  the  pseudo  Wigner  distribution  (with  L  =  Nw/2). 


The  concentration  measure  is  needed  in  the  algorithm,  Step  5. 
normalized  measure 


M\ 


where  Ad£  is  defined  in  [35]  as 


Here  we  use  the 


a<I=(EE  irKWY 

'  n  k  ' 


and  T(n,  k )  are  the  discrete  samples  of  a  non-negative  part  of  the  IRT. 


Example  1:  Consider  N  =  129  samples  of  noise  free  FM  signal  (4)  sampled  with 
t  =  nAt,  At  =  1/128,  n  —  0, 1, . . .  N  —  1.  The  signal  parameters  are  Ax  =  1,  fm  =  1.4, 
Arn  =  50  and  0rn  =  30°  The  spectrogram  and  the  S-method  of  the  considered  signal 
are  presented  in  Figure  1(a)  and  (b).  The  spectrogram  is  calculated  with  a  17-point 
Hann  window,  while  the  S-method  is  calculated  with  a  55  point  Hann  window  and 
L  =  8.  In  both  cases,  the  time-frequency  representation  is  calculated  at  each  available 
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time  instant  (i.e. ,  with  overlapping).  The  parameter  a  is  varied  from  0.2  to  15 
with  step  0.2.  For  each  a  ,  the  IRT  along  with  the  corresponding  concentration 
measure  is  calculated.  The  concentration  measure  /i(a)  is  presented  in  Figure 
1  (c)  and  (d),  where  the  minimum  of  this  measure  (i.e.,  highest  concentration)  is 
depicted  by  a  circle. 

Figure  1(e)  and  (f)  present  the  IRT  obtained  for  a  =  8.8.  The  maximum  position  in 
the  IRT  is  determined  and  the  modulation  parameters  are  estimated  in  the  spectro¬ 
gram  case  as 

frn  =  «  1-4,  Am  =  49.1,  em  =  31.4° 

Z7T 

and  in  the  S-method  case  as 

fm  =  «  1.4,  Am  =  48.1,  dm  =  32.2° 

Z7T 

As  we  can  see  in  both  cases,  the  modulation  parameters  are  very  close  to  true  values. 
It  means  that  the  method  will  not  be  too  sensitive  with  respect  to  the  time-frequency 
representation. 

The  IRT  for  a  =  7  (in  the  spectrogram  case)  and  a  =  10  (in  the  S-method  case)  are 
given  in  Figure  1(g)  and  (h).  These  subplots  illustrate  that  the  IRT  for  optimal  a 
(subplots  (c)  and  (d))  is  better  concentrated  than  the  IRT  with  another  a. 

Example  2:  The  estimation  procedure  is  applied  to  a  noisy  signal  s(nAt)  =  x(nAt)-\- 
e{n)1  where  the  noise  £{n )  is  a  complex  white  Gaussian  noise  with  SNR  =  0  dB.  The 
results  are  presented  in  Figure  2(a)  and  (b).  The  spectrogram  of  s(n )  is  presented 
in  Figure  2(a),  while  the  concentration  measure  of  the  IRT,  for  various  a,  is  given  in 
Figure  2(b).  In  the  spectrogram  calculation,  the  number  of  frequency  points  is 
101,  i.e.,  windowed  signal  is  zero-padded  prior  to  the  DFT  calculation.  Based  on 
the  IRT  obtained  for  optimal  a  (denoted  by  circle  in  Figure  2(b)),  the  modulation 
parameters  are  estimated  as 

fm  =  ^~  1-4,  Am  =  49.1,  em  =  30.2° 

Z7T 

and  the  resulting  sinusoidal  modulation  is  plotted  over  the  spectrogram  with  a  black 
line  Figure  2(a).  The  estimated  parameters  are  very  close  to  the  parameters  estimated 
for  noisy  free  case. 

A  multicomponent  signal  composed  from  a  sinusoidally  FM  component  (the  same  as 
in  Example  1),  a  linear  FM  component,  and  a  constant  frequency  component, 

s(t )  =  x(t )  +  0.6  exp  (j407r(f  —  0.8)2)  +  0.6  exp(j507rt) 
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Figure  1:  Modulation  parameters  estimation  for  the  mono-component  non-noisy  si¬ 
nusoidally  FM  signal.  Time  frequency  representation  (a),  (b);  concentration  measure 
(c),  (d);  inverse  Radon  transform  with  highest  concentration  (e),  (f)  and  inverse 
Radon  transform  when  the  parameter  a  is  not  optimally  chosen  (g),  (h). 
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Figure  2:  Parameter  estimation  of  the  noisy  signal  with  SNR  =  0  dB  (a),  (b)  and  the 
multicomponent  signal  (c),  (d).  Estimated  modulation  is  plotted  with  black  line  over 
the  spectrogram  image. 


is  considered.  The  results  obtained  with  the  proposed  procedure  are  presented  in 
Figure  2(c)  and  (d).  Estimated  parameters  are 

fm  =  «  1.4,  Am  =  49.1,  em  =  30.2°. 

2  7T 

From  this  example,  we  see  that  the  proposed  method  is  robust  to  the  noise  and  some 
other  possible  interferences. 

2.4  Multicomponent  signal  analysis 

This  approach  may  be  generalized  to  a  multicomponent  signal 

K  (  \ 

x(t)  =  Y  Ax]  exP  sin(27T  f^t  +  0{£)  )  +  e(t),  (5) 

k= 1  V  Jm  / 

where  e(t)  denotes  disturbing  components  and  noise.  Two  scenarios  are  possible.  One 
is  that  in  the  application  of  the  previous  algorithm,  in  Step  5,  the  concentration  meas¬ 
ure  p  of  the  obtained  IRT  produces  at  once  several  or  all  K  values  of  a  with  visible 
and  distinguishable  concentration  measure  peaks.  Then,  these  values  are  associated 
to  the  corresponding  signal  parameters,  as  in  Steps  6-9.  However,  due  to  different 
amplitudes  and  different  number  of  periods  in  the  time-frequency  plane  usually  only 
the  strongest  component  is  visible  in  the  concentration  measure.  In  this  case,  its 


DRDC-RDDC-201 4-R6 


9 


parameter  a  is  estimated  as  in  Step  5.  Other  parameters  are  estimated  for  this  com¬ 
ponent  as  in  Steps  6-9.  The  strongest  component  is  removed  and  the  algorithm  is 
used  on  the  remaining  signal  components,  until  the  energy  of  the  remaining  signal  is 
negligible.  After  parameters  of  all  components  are  found,  they  can  be  readjusted  by 
a  mean-squared  comparison  with  the  original  signal. 

Example  3:  Let  us  consider  a  multicomponent  noisy  signal  consisted  of  K  =  3 
sinusoidally  FM  components  of  the  form  (5).  Signal  parameters  are: 


a£1}  =  1,  /£>  =  1.4,  =  50,  =  30°, 

Ai2)  =  0.7,  /i2)  =  1,  A%!  =  35.7,  0%  =  -60°, 

A®  =  0.7,  /i3)  =  0.8,  At]  =  28.6  and  =  180°. 


The  proposed  method  estimates  parameters  of  one  component,  as  presented  in  Fig¬ 
ure  3(a)  and  (b).  The  estimated  parameters  are 


29.6° 


From  Figure  3(a),  we  can  see  that  estimated  modulation  parameters  highly  corres¬ 


pond  to  the  component  instantaneous  frequency.  Now  we  will  filter  out  the  estimated 
component . 


In  the  filtering  procedure,  the  original  signal  is  demodulated 


The  DFT  of  the  demodulated  signal  Xd(k)  =  DFT[xd(n)]  is  calculated  and  DC 
component  is  removed  by  putting  zero  value  to  Xd(0),  and  several  neighboring  points 
Xd(l),  Xd(N),  Xd(2),  Xd(N— 1),  etc.  Here  the  signal  length  is  N  =  129  and  we  remove 
7  points.  The  filtered  signal  is  obtained  by  the  inverse  DFT  Xf(n )  =  IDFT[Ad(/c)]. 

Finally  filtered  signal  is  modulated  in  order  to  cancel  frequency  shifts  in  the  remaining 
components  caused  by  the  demodulation 


Now  we  can  repeat  the  estimation  procedure  with  x(n)  =  xm(n )  and  estimate  the 
second  component  parameters.  The  results  are  presented  in  Figure  3(c)  and  (d).  The 
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Figure  3:  Multicomponent  signal.  Estimation  of  the  first  component  (a),  (b);  the 
second  component  (c),  (d);  and  the  third  component  (e)  and  (f).  Each  component  is 
removed  from  the  signal  after  the  estimation,  according  to  the  described  procedure, 
prior  to  next  component  estimation. 
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Figure  4:  Nonsinusoidal  modulation.  Triangulary  modulated  signal  (a),  (b);  Signal 
with  nonsinusoidal  modulation  and  varying  amplitude  (c),  (d). 

estimated  modulation  parameters  are 

fS}  =  I*  «  1-02,  ~  35.4,  9®  =  -63.4° 

In  the  next  step,  we  filter  the  estimated  component  and  proceed  to  the  parameters 
estimation  for  the  last  component.  Results  are  given  in  Figure  3(e)  and  (f)  and  the 
estimated  parameters  are 

/i3)  =  h  ~  °-796,  i-)  =  28-5’  ^)  =  -i78-7°- 

The  agreement  with  the  true  parameters  is  high. 


2.5  Nonsinusoidally  modulated  signals 

The  presented  estimation  procedure  could  be  used  even  if  the  analyzed  signal  is 
periodic,  but  not  sinusoidally  modulated.  We  will  illustrate  this  application  on  an 
example. 


Example  4:  Considered  a  triangularly  modulated  signal  x\  (t)  and  nonsinusoidal 
periodic  modulated  signal  X2(t)  with  a  varying  amplitude 


Xi(t)  =  exp  (  j  /  200  arcsin(cos(3.67ra))du 


x2(t)  =  A(t)  exp  (j  300\/arcsin(cos(3.67ra))(iu 
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where 


A(t)=exp(-(^)2). 


Although  the  proposed  method  is  derived  having  in  mind  the  sinusoidal  modulation, 
the  results  presented  in  Figure  4  clearly  show  that  the  applicability  of  the  proposed 
method  is  not  limited  to  the  sinusoidal  modulation  patterns  only.  The  estimated 
modulation  parameters  for  signal  X\  (t)  are 

fm  =  ^  ~  1-78,  Am  =  41.4,  6m  =  6.1° 

Z7T 

and  for  signal  X2  (t) 

L  =  ^  «  1-81,  Am  =  54.5,  9m  =  1.33°. 

Z  7T 

They  agree  with  fm  =  1.8  in  the  considered  signals.  The  closest  estimated  sinusoids 
are  presented  in  this  figure  as  well. 

2.6  Analysis  using  partial  data 

Assume  that  not  all  signal  samples  is  available.  In  this  case,  we  can  calculate  the 
spectrogram  only  at  time  instants/intervals  when  signal  samples  x{n)  are  available. 
This  procedure  will  be  illustrated  with  an  example. 

Example  5:  Consider  the  signal  defined  in  Example  3,  and  assume  that  samples  from 
intervals  20-28,  50-80,  and  95-110  are  missing.  Since  the  total  number  of  samples  is 
129  we  have  43%  of  missing  samples.  The  estimation  results  obtained  by  using  the 
available  signal  values  are  presented  in  Figure  5.  Regions  with  unavailable  samples  are 
presented  with  white  color  in  this  figure.  The  parameters  estimated  using  available 
samples  are 

«  1-401,  A%  =  49.4,  ^=30.8° 

«  0.986,  A$  =  35.0,  0 ^  =  -58.3° 

«  0.796,  A(A]  =  28.5,  =  -178.7°. 

Even  with  a  reduced  number  of  available  signal  values,  the  presented  method  pro¬ 
duced  accurate  estimates. 


fm 

J  m 

fm 

J  m 
7(3) 

J  m 


6.4 

27T 

6.4 
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Figure  5:  Multicomponent  sparse  signal.  First  component  estimation  (a),  (b);  Second 
component  (c),  (d);  Third  component  (e),  (f).  Missing  values  in  time-frequency 
representations  (a),  (c)  and  (e)  are  presented  in  white  color. 
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2.7  Experimental  data:  rotating  antenna  in  SAR 

Radar  returns  were  collected  from  a  rotating  antenna  using  a  APY-6  radar  in  a  SAR 
scenario  [2],  Using  these  data  sets,  the  m-D  features  relating  to  a  rotating  antenna 
were  extracted.  The  m-D  features  for  such  rotating  targets  may  be  seen  as  a  sinusoidal 
phase  modulation  of  the  SAR  azimuth  phase  history.  The  phase  modulation  may 
equivalently  be  seen  as  a  time-varying  Doppler  frequency. 

Figure  6a  shows  the  original  SAR  image  and  Figure  6b  displays  the  zoomed  in  SAR 
image  between  the  range  cells  115  and  130.  The  Doppler  smearing  due  to  the  rotating 
parts  is  often  well  localized  in  a  finite  number  of  range  cells.  It  is  reasonable  to  process 
the  Doppler  signal  for  each  range  cell  independently.  Since  the  prior  information 
about  the  location  of  the  target  is  known,  the  data  at  the  range  cells  123  and  124 
were  analyzed.  The  micro-Doppler  toolbox,  which  is  presented  in  Section  4  was  used 
to  estimate  the  rotation  rate  of  the  antenna.  Figure  7  illustrates  the  m-D  parameter 
estimation  results.  Figure  7a  shows  the  rotation  rate  at  range  cell  123  and  Figure  7b 
shows  the  rotation  rate  at  range  cell  124  The  estimated  rotation  rate  is  4.8  seconds, 
which  is  very  close  to  the  actual  value  of  4.7  seconds. 

2.8  Experimental  data:  hovering  helicopter 

The  experimental  data  used  in  the  analysis  that  follows  are  of  a  hovering  helicopter. 
For  a  helicopter,  the  main  rotor  blades,  the  tail  rotor  blades  and  the  hub  have  unique 
signatures  suitable  for  target  identification  [39-41],  Generally,  radar  returns  from  a 
helicopter  are  back-scattered  from  the  fuselage,  the  rotor  blades,  the  tail  blades  and 
the  hub  among  other  structures.  The  motion  of  the  rotor  blades  depends  upon  the 
coupling  between  the  aerodynamics  and  the  rotor  dynamics.  Each  blade  is  a  rotating 
aerofoil  having  bending,  flexing  and  twisting  motion.  The  radar  cross  section  of  a 
segment  of  the  blade  depends  upon  its  distance  from  the  centre  of  rotation,  its  angular 
position  and  the  aspect  angle  of  the  rotor  [40].  For  simplicity,  the  rotor  blade  can  be 
modelled  as  a  rigid,  linear  and  homogeneous  rod  rotating  about  a  hxed  axis  with  a 
constant  rotation  rate. 

The  rotational  motion  of  rotor  blades  in  a  helicopter  imparts  a  periodic  modulation 
on  radar  returns.  The  rotation-induced  Doppler  shifts  relative  to  the  Doppler  shift  of 
the  fuselage  (or  body)  occupy  unique  locations  in  the  frequency  domain.  Whenever 
a  blade  has  specular  reflection  such  as  at  the  advancing  or  receding  point  of  rotation, 
the  particular  blade  transmits  a  short  flash  to  the  radar  return.  The  rotation  rate  of 
the  rotor  is  directly  related  to  the  time  interval  between  these  flashes.  The  duration 
of  a  flash  is  determined  by  the  radar  wavelength  and  by  the  length  and  rotation  rate 
of  the  blades.  A  flash  resulting  from  a  blade  with  a  longer  length  and  a  radar  with  a 
shorter  wavelength  will  have  a  shorter  duration  [1], 
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Figure  6:  Top  -  The  original  SAR  image;  Bottom  -  The  zoomed  in  SAR  image. 


16 


DRDC-RDDC-201 4-R6 


Figure  7:  Rotation  rate  of  the  antenna  (a)  at  the  range  cell  123  and  (b)  at  the  range 
cell  124. 
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Figure  8:  (a)  Rotation  rate  of  the  main  rotor  blades  and  (b)  Rotation  rate  of  the  tail 
rotor  blades. 


18 


DRDC-RDDC-201 4-R6 


The  helicopter  employed  in  the  experiment  is  hovering  above  the  ground  at  a  height 
of  approximately  60  m  and  at  a  range  of  2.5  km  from  the  radar.  The  main  rotor 
comprised  five  blades  and  the  tail  rotor  consists  of  six  blades.  The  nr-D  toolbox  was 
used  to  estimate  the  rotation  rates  of  the  main  and  tail  rotor  blades  of  the  helicopter. 
The  rotation  rate  of  the  tail  rotor  blades  is  1030  rpm  (Figure  8a).  The  rotation  rate 
of  the  main  rotor  blades  is  203  rpm  for  this  helicopter  (Figure  8b).  Both  estimations 
are  in  agreement  with  the  actual  values. 

3  Micro-Doppler  removal  in  the  radar 
imaging  analysis 


The  nr-D  effect  appears  in  the  inverse  synthetic  aperture  radar  (ISAR)  imaging  when 
a  target  has  one  or  more  fast  moving  parts  [2, 5,  7-9].  Similar  effect  appears  in  the 
synthetic  aperture  radar  (SAR)  imaging,  as  well,  [2,49]. 

In  this  report,  we  use  the  L-statistics  based  approach  for  the  rigid  body  separation. 
In  order  to  remove  the  m-D  effect,  we  perform  the  TF  analysis  within  the  coherent 
integration  time  (CIT).  In  our  previous  approach  [5],  we  used  order  statistics  and 
several  TF  representations  with  various  windows.  The  obtained  TF  representations 
were  then  used  to  make  decision  whether  a  component  belongs  to  the  rigid  body  or  to 
the  fast  moving  target  point.  Here,  we  use  only  one  window  function  in  the  analysis. 
Order  statistics  is  performed  based  on  the  spectrogram,  while  the  rigid  body  signal 
synthesis  is  done  in  the  complex  TF  domain.  This  approach  is  very  simple  to  use 
and  produces  better  results  than  the  other  approaches.  It  is  also  robust  to  the  noise 
influence,  since  it  uses  the  L-statistics,  being  known  as  a  robust  signal  processing 
tool  [51].  The  L-statistics  application  to  the  complex  STFT  leads  to  a  form  of  super¬ 
resolution  representation,  as  well.  It  can  separate  very  close  rigid  body  components, 
even  when  that  is  not  possible  by  using  the  standard  Fourier  transform  (FT)  over  the 
entire  CIT.  The  proposed  method  can  be  easily  adapted  for  efficient  compensation  of 
a  residual,  uncompensated,  rigid  body  acceleration  in  the  presence  of  the  m-D  effects. 

In  order  to  improve  the  calculation  efficiency,  we  have  proposed  a  procedure  to  es¬ 
tablish  whether  there  is  any  target  return  in  a  considered  range  bin.  Moreover,  by 
bearing  in  mind  that  in  the  ISAR/SAR  analysis  only  some  range  bins  may  contain 
the  m-D  effect,  while  most  of  the  range  bins  are  m-D  free,  in  this  report,  we  have 
defined  a  criterion  for  detecting  ranges  which  contain  the  m-D  effects.  The  m-D 
removal  procedure  could  be  performed  only  for  these  particular  range  bins. 
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3.1  Time-frequency  analysis  and  L-statistics 


A  simpler  way  to  localize  the  signal  behavior  in  shorter  intervals,  within  the  CIT,  is 
in  applying  a  window  function  to  the  standard  FT.  The  resulting  short-time  Fourier 
transform  (STFT)  is  defined  as 


STFT(t ,  ft) 


( t)w(t 


t)e~jnrdT, 


(6) 


or  in  a  discrete  form 


M—l 

STFT(m ,  k)  =  s(i)w(i  -  m)e~j2nik/M,  (7) 

i= 0 

where  w(i)  is  a  window  function  used  to  truncate  the  considered  signal.  The  squared 
absolute  value  of  the  STFT  is  called  the  spectrogram.  In  most  of  the  provided 
examples,  we  will  use  a  Hanning  window.  The  window  width  is  Mw,  w{i)  ^  0  for 
—Mw/2  <  i  <  Mv,/ 2  —  1.  In  our  applications,  the  window  is  zero  padded  up  to 
M,  the  same  number  of  samples  as  in  FT,  so  that  we  have  the  same  frequency  grid 
in  the  STFT  as  in  the  FT.  Then,  we  can  later  easily  reconstruct  the  FT,  without 
interpolation,  with  the  concentration  close  or  equal  to  the  concentration  of  the  original 
FT.  We  know  that,  by  using  a  lag  window  w(i)  in  the  STFT,  the  concentration  in 
frequency  is  reduced,  as  compared  to  the  original  FT.  For  example,  if  the  lag  window 
width  is  Mw,  then  the  concentration  of  a  sinusoidal  signal  is  reduced  M/Mw  times, 
i.e.,  the  STFT-based  ISAR/SAR  image  of  a  rigid  body  point  (the  main  lobe  of  the 
FT  of  a  sinusoid)  would  be  approximately  M/Mw  times  wider  than  the  original  FT 
based  image  of  the  same  reflecting  point.  We  will  also  refer  to  this  effect  as:  the 
concentration  being  M/Mw  times  lower  in  the  STFT  than  in  the  original  FT. 

3.1.1  Restoring  the  high  FT  concentration  from  the  STFT 

The  concentration  could  be  restored  to  the  original  one  by  summing  all  the  low 
concentrated  STFT  (complex)  values  over  m.  Since  we  calculated  ST  FT(rn.  k)  with 
the  window  of  the  width  Mw .  there  are  two  possibilities  for  its  summation:  (a)  For 
all  time  instants  0  <  m  <  M  —  1,  when  the  signal  s{i)  has  to  be  zero-padded  for 
—Mw/ 2  <  i  <  0  and  M  <  i  <  M  +  Mwj 2  —  1;  (b)  For  instants  Mwj 2  <  m  < 
M  —  Mw/ 2,  when  zero-padding  of  s(i)  is  not  used.  Reconstruction  formula,  for  the 
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case  when  the  signal  is  not  zero-padded,  is 


M—Mw/2 

^  STFT{m,  k)  = 

m=Mw  /2 


M— 1 


M—Mw/2 

w(i  —  m) 

m=Mw/  2 


M—l 

=  s(i)Wl(i)e-j2,rik/M  = 
2=0 


^—j2,Kik/M 


SWl(k). 


(8) 


In  the  case  when  the  STFT  is  calculated  for  each  time  instant  (time  step  one  in 
the  STFT  calculation),  the  resulting  window  w\(i)  is  constant,  w\ (i)  =  const,  for 
Mw  —  1  <  i  <  M  —  My, ,  for  any  window.  It  means  that  during  the  most  of  the  CIT 
interval  we  have  the  normalized  resulting  window  W\{i)  being  close  to  the  rectangular 
one,  with  a  small  transition  at  the  ending  Mw  points.  The  FT  of  the  window  obtained 
during  the  process  of  reconstruction  produces  a  concentration  very  close  to  the  full 
range  rectangular  window  case  (i.e.,  no  window).  It  means  that  we  will  be  able 
to  reconstruct  the  FT  with  a  concentration  close  to  the  one  in  the  original  FT,  by 
using  low  concentrated  STFTs,  calculated  with  narrow  windows.  In  this  way,  we  will 
restore  the  high  concentrated  radar  image,  although  we  used  low  concentrated  STFT 
in  the  analysis.  The  transition  at  the  ending  points  of  w\  (z)  can  be  easily  overcome  by 
zero  padding  the  analyzed  signal  s{i)  with  Mw/2  samples  on  both  sides,  (as  explained 
before).  Then,  the  pure  rectangular  window  W\(i)  would  be  obtained,  for  any  window 
w(i).  The  analysis  is  not  restricted  to  the  step  one  in  the  STFT  calculation.  The 
same  resulting  window  would  be  obtained  for  a  step  equal  to  a  half  of  the  window 
width  [Mwj 2)  and  a  Hann(ing),  Hamming,  triangular  or  rectangular  window.  The 
same  is  valid  for  steps  equal  to  Mwj 4,  Mw/ 8,  etc. 

In  order  to  explain  how  this  mechanism  of  restoring  the  original  concentration,  by 
summing  low-concentrated  images,  works,  consider  s(t)  =  exjp(ju0t).  Its  FT  is  a 
delta  pulse  S'(fi)  =  2n <5(fi  —  ojq).  The  STFT  of  this  signal  produces 


STFT(t,  fi)  =  W(Q,  -  u o)  exp (j{n  -  u ;0)t).  (9) 


Let  us  now  analyze  the  result  of  summing  the  STFT  values  over  t: 

For  O  =  Co?o,  constant  values  of  W(0)  will  be  integrated  over  an  infinite  time  interval, 
with  the  phase  exp() (Q  —  oJo)t)  =  exp(jO),  producing  STFT (t,Q)dt  — >  oo  for 

fl  —  (jJ  o- 

For  any  other  H  not  equal  to  uq,  i.e.,  when  H  =  ujq  +  6,  9  ^  0,  we  will  have  the 
integration  W(6)  exp(jdt)dt  =  W{6)  exp(jdt)dt  =  0.  Therefore,  all  values 
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for  Vt  ^  ujq  are  averaged  out  to  zero  by  summation  of  the  low  concentrated  STFTs 
over  time. 

The  discrete  form  of  (9),  is 

M—l 

m  =  £  W(k  -  k0)ej2nm(k-k°yM  (10) 

m= 0 

with  the  same  conclusions  as  in  the  continuous  case.  Values  of  S(k),  when  the  signal 
is  not  zero-padded,  are  close  to  (10). 

3.1.2  Basic  idea  for  the  separation  of  a  rigid  body  and  fast 
rotating  part 

The  presented  mechanism  of  restoring  the  original  concentration  of  the  FT,  in  con¬ 
junction  with  the  knowledge  of  the  TF  behavior  patterns  of  fast  moving  and  rigid 
scattering  points,  lead  us  to  an  algorithm  for  the  m-D  free,  highly  concentrated,  radar 
image.  The  rigid  body  and  the  fast  moving  points  behave  differently  in  the  TF  rep¬ 
resentation  of  the  returned  radar  signal,  within  the  CIT.  The  rigid  body  signal  is 
almost  constant  in  time  (stationary),  while  the  fast- varying  m-D  part  of  the  signal  is 
highly  non-stationary.  This  part  of  signal  keeps  changing  its  position  in  the  frequency 
direction. 

For  the  illustration,  let  us  assume  that  the  signal  is  returned  from  one  point  of  a 
rigid  body  scatterer  and  one  point  of  a  fast  rotating  (m-D)  scatterer.  We  analyze  two 
cases  with  different  strengths  of  the  m-D  reflection.  In  the  first  case,  the  reflection 
coefficient  of  the  rigid  body  is  aB  =  1,  while  the  reflection  coefficient  of  the  fast 
moving  scatterer  is  aR  =  0.8.  The  STFT  representation  of  the  resulting  signal  is 
shown  in  Figure  9(a).  The  second  case  is  with  a  strong  m-D,  (Tr  =  15  and  the  same 
aB  as  in  the  previous  case.  The  STFT  representation  of  this  signal  is  shown  in  Figure 
9(e).  In  both  cases,  the  rigid  body  part  is  at  a  constant  frequency  for  all  t  within  the 
CIT,  while  the  fast  rotating  part  changes  frequency.  If  we  perform  sorting  over  the 
time  axis,  as  in  Figure  9(b,f),  we  will  not  change  the  result  of  the  summation  in  (8) 
since  it  is  a  commutative  operation.  By  summing  the  STFT  values  over  time,  from 
either  of  these  two  plots,  presented  in  Figure  9(a,b)  or  Figure  9(e,f),  we  will  get  the 
original  FT  of  the  corresponding  signal  Figure  9(c,g).  Note  that  any  value  of  (Tr  from 
(and  including  the  case  without  m-D)  aR  =  0  up  to  a r^>  a  B  will  not  significantly 
change  the  pattern. 

The  basic  idea  for  separating  the  rigid  body  and  the  fast  rotating  part  is  in  the  sorting 
of  STFT  values  of  the  returned  radar  signal  along  the  time  axis,  within  the  CIT.  Since 
the  rigid  body  return  is  stationary,  the  sorting  procedure  will  not  significantly  change 
the  distribution  of  its  values.  However,  the  fast-varying  m-D  part  of  the  signal  is 


22 


DRDC-RDDC-201 4-R6 


(a) 


frequency 


frequency 


frequency 


frequency 


frequency 


100  150  200  250 
frequency 


Figure  9:  Simulated  radar  signals  which  correspond  to  a  rigid  body  reflector  with 
<7b  =  1  and  a  rotating  reflector  with  reflection  coefficient  gr  =  0.8:  a)  Absolute  value 
of  the  STFT,  b)  Sorted  STFT  values,  c)  The  original  FT,  and  d)  The  reconstructed 
FT.  Simulated  radar  signals  which  correspond  to  a  rigid  body  reflector  with  gb  —  1 
and  a  rotating  reflector  with  reflection  coefficient  gr  =  15:  e)  Absolute  value  of  the 
STFT,  f)  Sorted  STFT  values,  g)  The  original  FT,  and  g)  The  reconstructed  FT. 


DRDC-RDDC-201 4-R6 


23 


highly  non-stationary,  occupying  different  frequency  bins  for  different  time  instants 
(in  the  case  of  flashes  it  exists  for  some  time-instants  only).  Its  existence  is  short  in 
time,  for  each  frequency,  over  a  wide  range  of  frequencies.  Thus,  after  sorting  the 
STFT  along  the  time  axis,  the  m-D  part  of  the  signal  has  strong  values  over  a  wide 
frequency  range,  but  for  a  few  samples  only.  By  removing  several  large  amplitude 
values  of  the  sorted  STFT,  for  each  frequency,  we  eliminate  most  or  all  of  the  m-D 
part  of  the  signal.  Summing  the  rest  of  the  STFT  values  over  time  we  will  get  the 
rigid  body  radar  image.  The  sinusoidal  m-D  pattern,  presented  in  Figure  9,  is  just 
an  example  of  such  a  signal.  This  idea  can  be  applied  on  any  non-stationary  signal 
form.  The  m-D  part  of  a  signal  is  non-stationary  by  definition. 

Let  us  consider  a  set  of  M  (or  M  —  Mw  if  the  signal  is  not  zero-padded)  elements  of 
the  STFT,  for  a  given  frequency  k, 

Sfc(ra)  =  {STFT(m,  k),m  =  0,1, ...,  M  —  1}. 

After  sorting  Sk(jn)  along  the  time,  for  a  given  frequency  k,  we  obtain  a  new  ordered 
set  of  elements  'Lfc(m)  G  Sk(m)  such  that  |\Et(0)|  <  |\Et(1)|  <  ....  <  |T/C(M  —  1)|.  Of 
course,  the  addition  is  commutative  operation,  so  if  we  use  the  whole  set,  we  get 

M— 1  M—l 

STFT(m,  k)=Y^  =  S(k)- 

m= 0  m= 0 


In  the  L-statistics  form  of  this  summation  we  will,  for  each  k ,  omit  M  —  Mq  of  the 
highest  values  of  Tfc(m)  and  produce  the  L-estimate  of  S(k),  denoted  by  Si{k),  as 

Mq-  1 

SL(k )  =  ^2  ^ k(m )  (11) 

m= 0 

where  Mq  =int[M(l  —  Q/100)]  and  Q  is  the  percent  of  omitted  values. 

To  illustrate  this  procedure,  we  eliminated  40%  of  the  top  amplitude  values  of  the 
STFT  from  the  previous  example.  In  this  way,  we  completely  eliminated  the  m-D 
component  from  the  TF  representation.  We  are  left  with  60%  of  the  low  amplitude 
values  of  the  STFT,  that  contain  only  the  rigid  body.  The  FT  reconstruction  is 
performed  based  on  these  values  only.  The  reconstructed  FTs  for  the  cases  of  a  weak 
and  a  strong  m-D  are  shown  in  Figure  9(d,h),  respectively.  The  FT  of  the  rigid 
body  is  in  both  cases  successfully  reconstructed  by  summing  60%  of  the  sorted  STFT 
samples,  remained  after  m-D  separation.  Note  that  the  result  is  not  significantly 
influenced  by  the  value  of  aR,  since  the  points  corresponding  to  the  m-D  signature 
are  removed,  meaning  that  their  values  are  almost  not  important. 

In  the  data  analysis,  this  approach,  based  on  elimination  of  a  part  of  data,  before 
analyzing  the  rest  of  the  data,  is  known  as  the  L-statistics  [51]. 
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Analysis  of  the  Missing  Values 

Since  we  have  eliminated  some  of  the  TF  representation  (TFR)  values,  we  will  analyze 
the  influence  of  incomplete  sum  in  (8).  This  is  the  same  theory  like  the  L-statistics 
theory  applied  to  the  noisy  or  non-noisy  data,  [51]. 

Assume  that  only  points  in  m  G  Dk  are  used  in  summation: 

sr.(k)  =  y.  STFT(m,k ),  (12) 

mG-Dfc 

where,  for  each  k,  Dk  is  a  subset  of  {0, 1,  2 , M  —  1}  with  Mq  elements. 

Within  the  framework  of  the  previous  analysis,  it  means  that  there  will  be  a  highly 
concentrated  component  S(k)  surrounded  by  several  low-concentrated  values 
Yhm^Dk  STFT(m,  k).  Note  that  the  amplitude  of  STFT(m,k )  is  M  times  lower 
than  the  amplitude  S(k),  since  S(k)  is  obtained  as  a  sum  of  M  values  of  the  STFT. 
In  general,  by  removing  let  say  (M  —  Mq)  values  in  m,  we  will  get  one  very  highly 
concentrated  pulse,  as  in  S(k),  and  [M  —  Mq)  values  of  low- concent  rated  components 
of  the  type  STFT(m,  k),  being  spread  around  the  peak  of  S(k)  and  summed  up  by 
different  random  phases.  Only  the  peak  value  is  summed  in  phase.  Consider: 

1.  Case  for  k  =  k0  corresponding  to  the  position  of  the  rigid  body  point:  At  this 
frequency,  all  terms  in  the  sum  (10)  are  the  same  and  equal  to  PF(O).  Thus,  the 
value  of  SL(k)  does  not  depend  on  the  positions  of  the  removed  samples.  Its  value  is 
SL(k0)  =  MqW{  0). 

2.  Case  for  k  =  l  +  ko,  where  l  %  0:  Removed  terms  in  (10)  are  of  the  form 
x;(m)  =  W (J)e:>2nml/M .  They  assume  values  from  the  set  T;  =  {W(£)ej27rmZ//M,  m  = 
0, 1, 2, ...,  M  —  1},  with  equal  probability,  for  a  given  l.  The  statistical  mean  of  these 
values  is  E{x.i(m)}  =  0  for  l  ^  0,  resulting  in  E{Sl(1  +  &o)}  =  0. 

The  resulting  statistical  mean  for  any  k  is 

E{SL(k)}  =  MQW(0)6(k-ko). 

The  higher  order  statistical  analysis  of  this  process  could  be  performed  in  detail,  but 
it  is  out  of  the  scope  of  this  report.  Here,  the  influence  of  the  number  of  missing  points 
to  the  concentration  of  the  reconstructed  FT  will  be  illustrated  by  an  example,  Figure 
10.  Here  we  consider  a  constant  frequency  signal,  without  m-D.  Its  FT  is  calculated 
and  presented  in  Figure  10(a).  Then  the  FT  is  reconstructed  based  on  25%,  50%  and 
75%  of  the  low  amplitude  values  of  the  STFT  for  each  k.  We  can  see  that  even  by 
taking  a  small  number  of  STFT  points,  we  still  keep  a  strong  peak,  since  it  is  summed 
in  phase,  Figure  10(b), (c),(d). 
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Figure  10:  The  FT  of  a  sinusoidal  signal:  (a)  Original ,  (b)  Reconstructed  by  summing 
75%  of  the  smallest  STFT  values,  for  each  k,  (c)  Reconstructed  by  summing  50%  of 
the  smallest  STFT  values,  (d)  Reconstructed  by  summing  25%  of  the  smallest  STFT 
values. 

Noise  Influence 

It  is  well  known  [51]  that  the  L-statistics  is  a  tool  for  robust  time- frequency  analysis. 
The  robustness  comes  from  the  fact  that  the  L-statistics  based  calculation  avoids 
highest  values,  which  are  the  most  influenced  by  noise.  Therefore,  we  may  expect 
that  by  using  the  L-statistics  we  will  not  degrade  the  radar  imaging  performance  in 
the  case  of  noise.  By  using  the  L-statistics  we  will  eliminate  a  part  of  the  signal,  that 
is  summed  in  phase  in  the  FT,  but  we  will  also  eliminate  the  signal  values  that  are 
mostly  corrupted  by  noise.  Thus,  with  the  elimination  of  the  m-D  we  will  improve 
the  overall  performance  in  the  noisy  signal  cases,  as  well.  In  the  case  of  impulse  noise, 
we  may  expect  significant  improvement,  even  in  the  case  without  m-D,  with  a  pure 
rigid  body.  The  effect  of  noise  is  statistically  analyzed  within  the  simulation  study  of 
this  report. 

High-Resolution  Property 

Consider  the  L-statistics  application  on  two  very  close  rigid  body  points: 

=  e-j2(yBi-AyBi)vBtuJo/c  _|_  e-j2(yBi+&.yBi)uBtu0/c 

with  a  very  small  A ysi  so  that  in  the  FT  of  the  signal  ,S'(T1),  calculated  over  the  entire 
CIT,  we  can  not  distinguish  these  components.  Since  the  resolution  in  the  Doppler 
direction  is  Rdopp  =  27 r/Tc  it  means  2A yBi  ~  2tt/Tc. 
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It  is  surprising,  but  if  we  use  low  concentrated  STFT  and  the  proposed  L-statistics, 
we  will  be  able  to  separate  these  components.  The  STFT  of  these  components  is 


STFT(t,Q)  = 

W(tt  +  (yBi  -  A yBi))e-*slte-XvBi-A''B* 
+W(Q  +  {yBi  +  A  yBi))e-jSHe-^VBi+AvBi)t 


with  normalized  frequency  0  for  2uqOJB/c  =  1. 

Note  that  the  STFTs  of  the  components  are  phase  shifted  for  A <p(t)  =  2AyBit.  Even 
for  small  2A yBi  ~  2i t/Tc  the  phase  shift  changes  are  of  the  order  A (p(t)  ~  2tt/Tc  x  Tc. 
It  means  that  it  could  easily  change,  during  the  CIT,  between  0  and  i r  or  even  more. 
Then,  there  will  be  time  instants  in  \STFT(t,  Q)  when  the  individual  STFTs  are 
summed  in  phase,  i.e. , 

when  \STFT(t,Q)\  =  |VF(Q  +  (yBi  -  A yBi))  +  PF(D  +  (yBi  +  Aym))\  .  Then  the  sig¬ 
nal  components  can  not  be  separated.  However,  there  will  be  also  the  instants  in  the 
STFT  when  the  components  are  with  opposite  phase, 

\STFT(t,Q)\  =  |VF(H  +  (yBi  -  A yBi))  -  VF(H  +  (yBi  +  AyBi))\,  so  that  the  signal 
components  are  clearly  separated.  We  can  see  that  in  the  first  case  the  values  of 
\STFT(t,  Q)  will  be  higher  than  in  the  other  case.  By  using  the  L-statistics  ap¬ 
proach  the  higher  values  will  be  eliminated,  while  the  lower  values,  that  are  well 
separated,  will  remain.  Thus,  we  may  achieve  high  signal  resolution  by  using  the  low 
concentrated  STFT  and  the  L-statistics,  even  in  the  case  when  the  separation  is  not 
possible  in  the  original  FT  over  the  whole  CIT. 

3.1 .3  Adaptive  percentage  of  missing  values 

There  are  two  possible  approaches  to  establish  the  threshold  for  the  elimination  of 
the  m-D: 

The  first  approach  is  to  assume  a  fixed  threshold  for  the  entire  ISAR/SAR  image: 
for  example,  removing  Q[%]  highest  values  for  each  frequency,  by  knowing  that  this 
will  not  disturb  significantly  the  obtained  image. 

A  more  sophisticated  approach  is  to  calculate  the  adaptive  threshold  for  each  range. 
The  adaptive  threshold  can  be  obtained  based  on  the  L-statistics.  To  use  the  L- 
statistics  approach,  we  will  sort  the  STFT  values  for  a  given  range  and  a  given 
frequency.  If  there  is  a  m-D,  then,  after  sorting  the  STFT  values,  there  is  a  region  of 
an  increase  of  the  sorted  STFT  values,  Figure  11  (upper  part).  Thus,  if  we  sum  the 
sorted  STFT  values  over  frequency,  we  will  get  a  function: 


(13) 
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Figure  11:  Sorted  STFT  values,  \Fi,(m)\,  with  a  colorbar  (upper).  The  values  of 
A(m),  obtained  by  summing  the  squared  values  of  the  sorted  STFT  along  the  fre¬ 
quency,  with  the  adaptive  threshold  Rl  obtained  based  on  the  average  value  of  10% 
of  its  smallest  values  -  thick  horizontal  line  (lower). 
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Now,  we  can  find  the  reference  level  Rl  based  on  the  mean  of  10%  of  the  low  amplitude 

M/ 10—1 

samples,  i.e. ,  based  on  M/10  values  of  the  sorted  ^(m),  Rl  =  Thr  10 

m= 0 

with  Thr  being  a  parameter,  usually  from  Thr  =  2  to  Thr  =  10.  For  example,  Thr  =  2 
means  that  we  will  use  all  the  values  in  the  sorted  1' /.-  (m)  whose  squared  values  are  up 
to,  for  example,  2  times  greater  than  the  mean  of  10%  lowest  squared  values.  Then  Q 
is  found  as  the  percent  of  A{m)  below  Rl-  In  the  case  when  there  is  no  m-D  it  means 
that  we  will  use  all  the  values,  since  the  stationary  values  are  close  to  the  mean  of 
the  lowest  values  for  all  time  instants.  If  there  is  a  m-D  then  the  value  outside  of 
the  stationary  points  will  start  to  increase  sharply  and  the  summation  will  stop.  The 
results  are  not  too  sensitive  to  these  values,  since  A(m )  is  a  fast  increasing  function 
when  the  m-D  starts  to  appear,  Figure  11  (lower  part). 

In  our  previous  work  [5],  we  have  used  the  L-order  statistics.  In  all  phases  of  the 
applications,  we  have  used  various  order  statistics  of  the  absolute  values  of  the  STFT 
with  various  window  widths.  Here,  we  use  only  one  window  function  for  the  analysis. 
Then,  after  sorting  the  absolute  STFT  values  and  defining  the  threshold  Q  (by  simply 
assuming,  for  example  Q  =  50%,  or  calculating  an  adaptive  threshold),  we  return  back 
to  perform  all  the  calculations  in  the  complex  STFT  domain.  In  this  way,  we  obtain 
a  very  simple  and  efficient  model  for  the  calculation,  while  the  results  are  improved 
with  respect  to  those  obtained  by  the  procedure  proposed  in  our  previous  work. 

3.1 .4  Algorithm  for  the  micro-Doppler  effects  removal 

The  simplest  way  to  use  the  proposed  method  is  in  applying  the  L-statistics  approach 
to  all  range  bins,  with  a  constant  threshold,  for  example  Q  =  50%.  In  this  case  the 
m-D  will  be  separated,  while  the  rigid  body  will  not  be  degraded.  In  the  case  of 
impulse  noise  we  will  benefit  from  this  procedure  in  each  bin.  Also,  if  there  are  close 
rigid  body  components,  then  this  approach  will  help  to  resolve  them,  as  it  is  described 
earlier.  However,  this  way  of  calculation  increases  the  calculation  complexity. 

In  radar  data  analysis,  only  a  part  of  the  range  bins  may  contain  the  m-D,  while  most 
of  them  will  be  m-D  free.  We  can  improve  the  calculation  complexity  by  defining  a 
procedure  to  avoid  processing  of  the  range  bins  where  we  can  conclude  that  there  is 
no  target  return  or  there  is  no  m-D  part  of  signal.  If  there  is  an  m-D  point,  then 
a  constant  or  adaptive  threshold  for  the  L-statistics  application  in  the  m-D  removal 
could  be  used  in  order  to  separate  the  m-D,  while  most  of  the  rigid  body  part  of 
signal  is  preserved.  If  there  is  no  m-D  in  the  considered  range  bin,  then  processing 
for  this  range  is  not  necessary.  If  the  FT  is  already  well  focused,  or  if  there  is  no 
returned  target  signal  in  the  considered  range  bin,  then  the  FT  can  be  used  as  it  is. 
Of  course,  this  classification  procedure  is  optional,  and  we  will  not  lose  anything  in 
the  radar  image  quality  if  we  apply  the  presented  method  in  some  bins  with  already 
focussed  image  or  where  there  is  no  target  return.  Thus,  the  thresholds  in  the  next 
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procedure  may  be  chosen  in  quite  a  conservative  way,  to  be  sure  not  to  miss  any  bin 
with  m-D,  allowing  false  m-D  detections. 


The  range  bins  classification  procedure  will  be  presented  as  an  algorithm. 
The  algorithm  consists  of  the  following  steps: 


Step  1  :  Detect  whether  there  is  a  returned  target  signal  in  the  considered  range  bin, 
by  using 

max{|S'(/c)|}  >  e,  (14) 

where  £  =  0.02max{|S'(fc,  l)\}  in  the  noiseless  case  or 

e  =  max{0.02  max{|S'(A:,  Z) | },  2y/cr/M}  in  the  presence  of  noise  (here  S(k,  l )  is  the  2D 
FT  of  the  received  signal,  i.e. ,  full  radar  image  matrix  for  all  ranges  and  cross- ranges, 
while  S(k)  is  the  radar  image  for  a  given  range  bin).  The  standard  deviation  o  of 
noise  in  the  radar  image,  for  a  given  range  bin,  can  be  estimated  as  [53] : 


medianj |Re{S'(/c)}  —  R e{S(k  —  1) } | ,  k  =  2, ..,  M } 

0.6745^ 


(15) 


for  the  real  part  of  S(k).  The  same  applies  for  the  imaginary  part.  If  there  is  no 
target  in  the  returned  signal,  do  not  perform  the  following  steps  and  take  S(k)  as  the 
radar  image  for  the  considered  range  bin.  If  there  is  a  target  signal,  continue. 


Step  2  :  Detect  whether  there  are  any  m-D  effects  in  the  considered  range  bin: 

max{|5'(fc)|} 

mean{|5'(A:)|} 

The  previous  relation  is  the  simplest  concentration  measure.  In  the  case  of  a  highly 
concentrated  signal,  R  is  close  to  M,  while  for  a  low  concentrated  signal,  this  value 
is  close  to  1.  We  say  that  the  signal  is  well  concentrated  (rigid  body  only)  if  R  >  10. 
This  threshold  was  successfully  tested  for  various  scenarios  and  multi-component 
signals.  If  the  signal  is  well  concentrated  take  its  FT  S(k)  as  the  radar  image  for  the 
considered  range  bin.  If  the  signal  is  not  well  concentrated,  continue. 


Step  3  :  Remove  the  m-D  effects  and  reconstruct  the  FT  of  the  reflectors  that  corres¬ 
pond  to  the  rigid  body  in  the  considered  range  bin.  Remove  the  values  of  the  STFT 
that  are  higher  than  the  threshold  and  for  each  frequency  sum  the  rest  along  the 
time.  Take  the  obtained  FT  S^fk )  (12)  as  the  radar  image  for  the  considered  range 
bin. 
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Figure  12:  (a)  The  STFT  of  a  signal  consisting  of  one  rigid  body  component  and  four 
sinusoidally  modulated  components,  (b)  The  sorted  STFT  of  the  same  signal,  (c) 
The  original  FT  of  the  signal,  (d)  The  FT  of  rigid  body,  reconstructed  by  summing 
the  STFT  values  remaining  after  sorting  and  eliminating  samples  that  correspond  to 
the  m-D  effect. 

3.2  Results 

Example  1:  The  proposed  method  is  tested  on  a  signal  with  one  rigid  body  point  and 
four  sinusoidally  modulated  components  (used  to  model  rotating  reflectors), 

K 

sign)  =  (Tb^2  exp  {jyBim}  (16) 

i= 1 
p 

+aR  ^2  exP  {j[ymim  +  Am  sin {umm  +  <^)]}  , 

1=1 

with  K  =  1,  P  =  4,  aB  =  1,  crR  =  3,  yBl  =  0.47T,  ARi  =  [96,  48,  64,  24],  uRi  =  7t/128, 
yRoi  =  7 r  and  ipi  =  0,  for  i  =  1,2,  3, 4.  The  STFT  of  this  signal  is  presented  in  Figure 
12(a).  The  m-D,  although  moderate,  significantly  covers  the  rigid  body,  i.e. ,  the  part 
of  the  constant  frequency  component  is  almost  invisible  in  the  sinusoidal  patterns. 
The  sorted  STFT  is  shown  in  Figure  12(b).  Then,  the  highest  STFT  values  are 
removed,  for  each  frequency  in  the  reconstruction  phase.  A  constant  threshold,  with 
Q  =  60%  is  used  here.  The  FT  reconstructed  from  the  remaining  STFT  samples  is 
shown  in  Figure  12(d).  The  rigid  body  is  successfully  reconstructed  in  the  presence 
of  the  m-D.  The  original  FT  of  the  analyzed  signal  is  given  in  Figure  12(c). 

Example  2:  Here,  we  analyze  a  signal  with  10  components:  K  =  5  components  with 
constant  frequency  (used  to  model  rigid  body  reflectors)  and  P  =  5  sinusoidally 
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Figure  13:  (a)  The  STFT  of  a  signal  consisting  of  five  rigid  body  components  and  five 
sinusoidally  modulated  m-D  components  (b)  The  sorted  STFT  of  the  same  signal,  (c) 
The  original  FT  of  the  signal,  (d)  The  FT  of  rigid  body,  reconstructed  by  summing 
the  STFT  values  remaining  after  sorting  and  eliminating  samples  that  correspond  to 
the  m-D  effect. 

modulated  components  (used  to  model  rotating  reflectors),  (16)  with:  aB  =  1,  crR  = 
15,  ym  =  [1.9tt,  1.95? r,  2vr,  2.05vr,  2.1? r],  Am  =  [150,  300,  200,  440,  200],  ujRi  =  [tt/256, 
7t/512,  7t/256,  7t/512,  7t/256],  yRoi  —  0  and  =  [0,  — 7r/ 3,  7r/6,  — 27t/3,  0],  for 
i  =  1,2,  3, 4,  5,  M  =  1024  and  Mw  =  64.  The  STFT  of  this  signal  is  shown  in  Figure 
13(a).  The  constant  components,  that  correspond  to  the  rigid  body,  are  not  well 
separated  in  the  TF  plane.  Moreover,  they  are  covered  by  the  sinusoidally  modulated 
patterns  which  represent  the  m-D  effects  of  the  rotating  reflectors.  If  we  sort  the 
STFT  values  along  time  axis,  then  the  representation  of  the  rigid  body  parts  does 
not  change,  since  it  is  constant  during  the  whole  CIT,  Figure  13(b).  On  the  other 
hand,  the  fast  rotating  parts  occupy  only  a  small  time  intervals  over  a  wide  region  of 
frequencies.  They  lie  in  high  value  regions  of  the  sorted  transform.  Thus,  they  will 
be  eliminated  by  removing  the  highest  STFT  values,  for  each  frequency.  An  adaptive 
threshold,  with  Thr  =  5  is  used  here.  The  results  are  not  sensitive  to  the  value  of 
Thr.  The  reconstructed  FT,  obtained  by  summing  the  rest  of  the  STFT  (12)  along 
the  time  is  shown  in  Figure  13(d).  We  can  clearly  see  five  peaks  that  correspond  to 
the  five  rigid  body  reflectors.  The  original  FT  is  shown  in  Figure  13(c).  It  cannot  be 
used  even  to  determine  the  number  of  components  in  the  analyzed  signal. 


Helicopter  Data  Analysis 

Example  3:  In  this  example,  we  first  present  a  new  simulation  approach  to  the  data 
of  a  German  Air  Force  Bell  UH-1D  Helicopter  known  also  as  ‘Iroquois’  presented  in  [5]. 
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Here,  the  simulation  is  performed  according  to  the  variable  flashing  reflection 
coefficients,  rather  than  just  by  using  a  mathematical  form  that  would  produce 
the  data  as  in  [6].  Several  effects  are  emphasized  in  the  TFR  Figure  14(a).  The 
stationary  patterns  along  the  time-axis  correspond  to  the  rigid  body  reflection. 
The  motion  of  two  main  blades  is  modeled  by  two  rotating  reflectors,  producing 
sinusoidal  FM  signals  with  a  large  magnitude  in  the  frequency  direction,  (18).  The 
main  rotor  flashes  are  simulated  by  signals  producing  lines  that  connects  extreme 
points  of  the  sinusoidal  FM  signal,  along  the  time  axis.  The  smaller  pulses  that  can 
be  seen  on  the  right-hand  side  of  Figure  14(a)  correspond  to  the  tail  rotor  flashes, 
and  they  are  simulated  here  by  taking  into  account  the  physical  meaning  of  its 
appearance.  Namely,  these  flashes  correspond  to  the  periodic  alignment  of  the  main 
and  tail  rotors  to  maximally  reflect  the  radar  signal  when  they  are  normal  to  the 
line-of-sight.  Therefore,  we  use  here  an  angle  dependent  reflection  coefficient 


a(t)  =  exp (—30  \sm(2nt/TRot)\),  (17) 

where  the  reflection  takes  value  1  when  t  =  kTRot/ 2  and  \sm(2irt/TRot)\  =  0,  while  for 
other  t,  30  \sm(2nt /TRot)\  assumes  high  values  and  the  reflection  coefficient  is  small. 
Note  that  other  effects  that  can  be  observed  in  a  radar  image,  including  multi-path, 
are  not  considered  here. 


The  simplified  model  of  the  reflected  UH-1D  signal  can  now  be  written  as 

s(t)  =  XRIG{t)  +  XROT{t )  +  XpL_M{t )  +  XFL_r{t), 

where  xRio(t),  xrot (£),  Ffl_m(£)  and  xfl_t( t)  represent  signals  caused  by  the  rigid 
body,  rotation  of  the  main  rotor,  and  the  main  and  tail  rotor  flashes,  respectively. 
The  signal  is  considered  within  the  interval  of  400  ms,  sampled  with  a  rate  of  At  = 
1/48  ms.  Four  sinusoidal  components,  caused  by  the  rigid  body,  are  at  the  frequencies 
—  10.3  kHz,  —2.5  kHz,  2.3  kHz  and  2.7kHz.  Two  components  at  —0.4kHz  and  0.4kHz 
correspond  to  the  modulated  time  tones  commonly  added  to  the  data  tape  [56] .  The 
sinusoidal  FM  signals,  corresponding  to  the  rotation  of  the  main  rotor  blades,  are 
modeled  as 


xROT{t)  =  aROT[ej2nAROTSin{2^t/TROT)  (18) 

^_e~j2TrAROT  sin(27rt/THoT)j 

where  cr rot  =  10,  Trot  =  175  ms  and  ARqt  =  529.19.  The  main  and  tail  rotor 
flashes  are  modeled  as 


128 

XFL_M(t )  =  2.5  ^ 

k= 1 


+  64  -30|sin(27rt/175)| 
128 


x  cos (25. 98 k  sin ( 


27 rt 
175 


)), 
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Figure  14:  (a)  The  STFT  of  a  simulated  signal  of  a  German  Air  Force  Bell  UH-1D 
Helicopter,  (b)  The  sorted  STFT  of  this  signal,  (c)  Original  FT  of  the  signal,  (d)  The 
FT  of  the  rigid  body,  reconstructed  by  summing  the  lowest  absolute  STFT  values, 
(e)  The  FT  of  the  rigid  body,  reconstructed  by  the  proposed  method. 


and 


128 

XFL_r{t )  =  2.5  ^2  e  >'e 

k= 64 


— 30|sin(27rt/35. 8)|  (j  (2.66/c  sin(47rt/35. 8))) 


The  signal  is  corrupted  by  a  moderate  Gaussian  noise.  To  compare  our  simulation 
with  the  real  one  (for  the  m-D  and  rigid  part  values  ratio)  refer  to  [6,56]. 


The  proposed  algorithm  for  the  rigid  body  separation  is  applied  to  the  simulated 
helicopter  signal.  The  sorted  STFT  is  shown  in  Figure  14(b).  We  can  see  that  the 
STFT  values  corresponding  to  the  rotating  parts  are  in  the  high  value  region.  The 
reconstructed  FT  is  shown  in  Figure  14(e).  All  5  reflectors  that  correspond  to  the 
rigid  body  are  successfully  recovered.  The  original  FT  is  presented  in  Figure  14(c), 
while  the  reconstructed  FT  obtained  by  summing  the  absolute  values  of  the  remaining 
STFT  samples  is  presented  in  Figure  14(d). 
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Figure  15:  (a)  The  STFT  of  a  signal  consisting  of  two  very  close  components  with 
constant  frequency  and  one  sinusoidally  modulated  component,  (b)  The  sorted  STFT 
of  the  same  signal,  (c)  The  original  FT  of  the  signal,  (d)  The  FT  of  rigid  body, 
reconstructed  by  summing  the  lowest  STFT  values. 

High-Resolution  Analysis 

Example  4:  Two  very  close  rigid  body  reflectors  in  the  presence  of  m-D  effects  are 
simulated  in  this  example 

s(m\  —  _|_  e~j205nm/M  _|_  IQgi58  cos(2™/M) 

where  M  =  256  samples  are  used.  The  window  with  Mw  =  32  zero-padded  to  M 
is  used  for  the  STFT  calculation.  The  STFT  of  the  analyzed  signal  is  presented  on 
Figure  15(a),  while  the  sorted  STFT  is  presented  on  Figure  15(b).  It  can  be  seen  in 
Figure  15(a),  that  there  are  time  instants  when  the  STFTs  of  the  close  components 
are  summed  with  opposite  phase,  and  they  appear  as  separated.  On  the  other  hand, 
when  the  close  components  are  summed  in  phase,  they  are  not  separated.  Moreover, 
as  it  can  be  seen  from  the  sorted  STFT,  presented  in  Figure  15(b),  when  the  STFTs 
are  summed  in  phase,  the  resulting  STFT  is  higher.  Consequently,  by  removing  the 
highest  values  of  the  STFT,  the  remaining  lower  values  are  well  separated;  thus,  the 
close  components  are  separated.  The  FT  reconstructed  by  summing  over  time  50%  of 
lowest  samples  of  the  STFT,  is  shown  in  Figure  15(d),  while  the  original  FT  is  shown 
in  Figure  15(c).  We  can  see  that  the  separation  of  the  close  components  is  achieved 
by  the  proposed  method,  although  it  is  not  possible  in  the  original  FT  (the  distance 
between  two  maxima  positions  is  biased). 
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Figure  16:  Mean  absolute  error  as  a  function  of  noise  variance  calculated  for  the  case 
of  one  rigid  body  reflector  without  m-D.  The  solid  line  corresponds  to  the  proposed 
method,  while  the  dashed  line  corresponds  to  the  full  FT. 

Noise  Influence  Analysis 

Example  5:  One  stationary  reflector  and  one  m-D  reflector  are  considered.  Complex 
valued,  white  Gaussian  noise  e(t),  with  variance  a2,  is  added 

S(m)  =  e-j°-75™  +  CTRei58coS(2™/256)  + 

where  (Tr  is  the  reflection  coefficient  of  the  m-D  reflector.  The  noise  variance  is 
varied  within  a  wide  range  0  <  a2  <  72  (from  the  case  without  noise  up  to  the 
case  when  noise  dominates),  with  step  1.  For  each  variance  value  from  this  range, 
1000  Monte  Carlo  simulations  are  performed.  In  each  realization,  we  have  found  a 
position  of  the  maximum  in  the  L-statistics  based  estimate  of  the  FT,  Sb(/c).  Then, 
the  error  is  calculated  as  a  difference  of  this  position  and  the  true  signal  frequency. 
The  mean  absolute  error  is  calculated  for  each  variance  for  1000  realizations  and  the 
mean  absolute  error  is  plotted  for  various  noise  variance  values.  For  the  rigid  body 
FT  reconstruction  we  used,  for  each  frequency,  50%  of  the  smallest  STFT  values  in 
the  L-statistics  summation. 

We  start  with  the  case  of  pure  stationary  point  ctr  =  0,  to  see  how  the  L-statistics 
approach,  with  50%  of  values,  influences  the  results.  It  is  well  known  that  the  FT 
transform  is  theoretically  the  best  (ML)  estimator  for  a  pure  sinusoid  in  Gaussian 
noise.  The  corresponding  mean  absolute  error  is  depicted  in  Figure  16.  The  solid  line 
corresponds  to  the  proposed  method,  while  the  dashed  line  corresponds  to  the  full 
FT.  The  FT  is  well  reconstructed  with  the  proposed  L-statistics  based  method  and 
the  estimation  results  are  not  degraded  with  respect  to  the  full  FT,  in  this  simple 
case,  when  the  FT  is  the  ML  estimator. 

We  have  analyzed  noise  influences  in  the  case  of  ctr  =  5,  as  well.  For  the  noiseless  case, 
the  STFT  is  shown  in  Figure  17(a),  while  the  sorted  STFT  is  shown  in  Figure  17(b). 
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The  original  FT  is  presented  in  Figure  17(c).  The  FT  reconstructed  by  summing, 
for  each  frequency,  50%  of  the  lowest  STFT  samples  is  shown  in  Figure  17(d).  The 
same  plots  for  the  case  of  at  =  4.5,  SNR  =  -6.53  dB  are  shown  in  Figure  17(e-h). 
The  signal  to  noise  ratio  (SNR)  is  calculated  as  the  rigid  body  part  of  the  signal 
to  the  noise  ratio,  in  all  cases.  We  can  see  that  the  proposed  method  successfully 
reconstructs  the  FT  of  the  rigid  body  in  the  presence  of  m-D  and  noise. 

The  mean  absolute  error  of  the  proposed  method  and  the  original  FT  is  shown  in 
Figure  17(i).  The  proposed  method  does  not  only  reconstruct  the  FT  successfully, 
but  also  eliminates  the  m-D  effect  and  outperforms  the  original  FT,  whose  estimation 
performance  is  degraded  by  the  m-D  effect. 

Example  6:  We  analyzed  one  more  case  of  one  stationary  reflector  and  one  m-D 
reflector  in  the  presence  of  noise.  Here,  the  m-D  reflector  is  stronger  and  closer  to 
the  stationary  one.  The  corresponding  signal  is  of  the  same  form  as  in  the  previous 
example,  but  with  =  10,  while  the  rigid  body  reflector  signal  component  is  at  the 
frequency  fs  =  0.125  Hz.  The  same  statistical  analysis  and  reconstruction  procedure 
as  in  the  previous  example  are  performed. 

For  the  noiseless  case,  the  STFT  is  shown  in  Figure  18(a),  while  the  sorted  STFT 
is  shown  in  Figure  18(b).  The  original  FT  is  presented  in  Figure  18(c).  The  FT 
reconstructed  by  summing,  for  each  frequency,  50%  of  the  lowest  STFT  samples  is 
shown  in  Figure  18(d).  The  same  plots  for  the  case  of  a2  =  4.5,  SNR  =  — 6.53dB 
are  shown  in  Figure  18(e-h).  We  can  see  from  Figure  18(d)  and  Figure  18(h)  that 
the  performance  of  the  proposed  method  does  not  degrade  even  in  the  case  of  strong 
m-D  reflector  positioned  close  to  the  rigid  body  reflector;  in  this  case,  the  stationary 
and  m-D  components  are  crossing  in  the  STFT,  Figure  18(a)  and  Figure  18(e).  The 
proposed  method  continues  to  successfully  reconstruct  the  FT  of  the  rigid  body  in 
the  presence  of  noise,  Figure  18(h),  while  the  FT  is  not  even  able  to  indicate  that 
there  is  a  rigid  body  reflector,  at  all,  Figure  18(g). 

The  mean  absolute  error  of  the  proposed  method  and  the  original  FT  are  shown 
in  Figure  18(i).  From  the  presented  statistics,  we  can  confirm  that,  even  in  the 
presence  of  noise  and  close  reflectors  with  strong  m-D  effects,  the  proposed  method 
successfully  reconstructs  the  FT  of  the  rigid  body,  while  the  original  FT  completely 
fails  to  indicate  the  rigid  body  existence. 

Non- Compensated  Rigid  Body  Acceleration 

Example  7:  In  this  case  an  accelerating  rigid  body  target  is  considered  and  examined. 
The  received  radar  signal  that  corresponds  to  an  accelerating  target  in  the  ISAR 
systems  is  a  linear  FM  signal.  Similarly,  in  SAR  systems  the  target  motion  may 
induce  linear  frequency  modulation  in  the  received  radar  signal  [9].  Therefore,  we 
simulated  three  rigid  body  reflectors  as  three  linear  FM  components  with  the  chirp- 
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Figure  17:  One  rigid  body  reflector  and  one  m-D  reflector.  Noiseless  case  (a -d):  (a) 
The  STFT  absolute  value,  (b)  Sorted  STFT,  (c)  The  original  FT,  (d)  Reconstructed 
FT.  A  realization  of  noisy  case  with  ni  =  4.5  (e  h):  (e)  The  STFT  absolute  value, 
(f)  Sorted  STFT,  (g)  The  original  FT,  (h)  Reconstructed  FT.  (i)  Mean  absolute  error 
as  a  function  of  noise  variance  in  1000  noisy  realizations. 


38 


DRDC-RDDC-201 4-R6 


<u 

13 

3 

a 

3 


<D 

13 

3 


s 

index 

frequency 

frequency 

o' 

amplitude 

0  50  100  150  200  250  0  50  100  150  200  250 

s 

'o' 

index 

frequency 

frequency 

SAVwJ 

_ * _ * _ vw% 

amplitude 

0  50  100  150  200  250 


0  50  100  150  200  250 


Noise  variance 


Figure  18:  One  rigid  body  reflector  and  one  close  m-D  reflector.  Noiseless  case 
(a -d):  (a)  The  STFT  absolute  value,  (b)  Sorted  STFT,  (c)  The  original  FT,  (d) 
Reconstructed  FT.  A  realization  of  noisy  case  with  rr2  =  4.5  (e-h):  (e)  The  STFT 
absolute  value,  (f)  Sorted  STFT,  (g)  The  original  FT,  (h)  Reconstructed  FT.  (i) 
Mean  absolute  error,  in  1000  noisy  realizations,  as  a  function  of  noise  variance. 
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rate  a.  In  order  to  show  that  our  algorithm  will  not  remove  only  the  m-D  induced 
by  vibrating  and  rotating  targets,  here  we  have  also  used  a  more  complex  form  of 
the  m-D.  In  this  example,  we  will  show  that  the  algorithm  is  robust  to  the  effects  of 
acceleration. 


The  STFT  of  the  analyzed  signal  is  presented  in  Figure  19(a).  We  can  clearly  see 
that,  as  a  result  of  the  acceleration,  the  TFR  of  the  rigid  body  part  of  the  signal  is 
not  stationary  during  the  time.  Consequently,  it  is  difficult  to  separate  it  from  the 
m-D  in  the  sorted  STFT,  Figure  19(b).  Namely,  if  we  perform  the  m-D  separation  by 
removing  50%  of  the  highest  STFT  samples,  as  we  did  in  the  examples  where  there 
was  no  need  for  the  motion  compensation,  we  would  reconstruct  the  FT  of  the  rigid 
body  as  presented  in  Figure  19(d).  Here,  we  have  removed  a  significant  part  of  the 
rigid  body  points,  as  well. 


In  the  analysis  of  the  rigid  body  with  uncompensated  acceleration,  we  should  first 
compensate  the  remaining  acceleration.  This  is  not  possible  in  the  original  signal, 
since  the  m-D  signatures  prevent  us  from  properly  compensating  the  remaining  ac¬ 
celeration.  However,  the  application  of  the  proposed  method  for  the  m-D  removal 
can  solve  this  problem,  as  well.  We  will  use  the  Local  Polynomial  Fourier  Transform 
(LPFT) 

poo 

LPFT(t ,  ft)  =  /  s(t)w(t  -  t)e-j^T+aT^dr,  (19) 


instead  of  the  STFT,  where  the  term  exp^—joiT2)  is  used  to  compensate  the  linear  fre  - 
quency  modulation  of  the  rigid  body  part  ofthe  signal,  LPFT(t,  ft)  =  FT{s{r)ew(T — 
f)}.  The  parameter  a  is  not  known  in  advance,  but  we  know  that  it  can  take  val 
ues  from  a  set  A  =  [— amax,  amax],  where  amax  is  the  chirp  rate  corresponding  to  the 
maximal  expected  acceleration  (positive  or  negative),  [57].  In  this  example  we  used 
A  =  [—2  :  0.25  :  2],  Now,  a  can  be  estimated  as  the  value  from  the  set  A  for  which  we 
obtain  the  highest  concentration  of  the  reconstructed  rigid  body  (compensated  FT) 
based  on  the  LPFT  and  the  L-statistics,  with,  for  example,  Q  =  50%.  The  recon¬ 
structed  FT,  by  using  50%  of  the  lower  LPFT  values,  will  be  denoted  by  SL,a(k)  . 
Its  concentration  is  calculated  using  the  time- frequency  concentration  measure  [58], 


H(a) 


M—l 

£  \suk)P 


(20) 


with  p  =  1.  The  LPFT,  calculated  with  the  estimated  optimal  value  of  a  =  1.25, 
which  results  from  H(a),  is  shown  in  Figure  19(e).  The  linear  frequency  modulation 
is  compensated  by  a  in  (19).  Thus,  with  optimal  a  we  have  components  with  almost 
constant  frequency  in  the  TFR  representation  of  the  rigid  body  reflectors.  In  this 
way,  we  have  successfully  reconstructed  the  rigid  body  and  removed  the  m-D  part, 
as  it  is  presented  in  Figure  19(h).  The  procedure  is  not  too  sensitive  to  a.  Very  good 
results  are  obtained  with  neighboring  values  a  =  1.0  and  a  =  1.5. 
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Note,  that  it  would  be  impossible  to  estimate  the  chirp-rate  a  from  the  original  signal, 
without  employing  the  proposed  algorithm  for  the  m-D  removal. 


Real  Data  Application 


frequency  frequency 


frequency 


0  50  100  150  200  250 


frequency 


Figure  19:  Accelerating  rigid  body  with  a  complex  form  of  the  mD.  (a)  TFR  of  the 
signal  without  motion  compensation,  (b)  Sorted  TFR  of  the  original  signal,  (c)  Ori¬ 
ginal  FT  of  the  analysed  signal,  (d )  Reconstructed  FT  of  the  accelerating  rigid  body 
without  motion  compensation,  (e)  TFR  of  the  signal  after  motion  compensation,  (f) 
Sorted  TFR  of  acceleration  compensated  signal,  (g)  The  FT  of  the  original  signal 
with  motion  compensation,  (h)  Reconstructed  FT  of  the  accelerating  rigid  body  with 
motion  compensation. 

Example  8:  The  proposed  algorithm  is  tested  on  real  data  in  this  example.  The 
examined  data  were  collected  using  an  X-band  radar  operating  at  9.2  GHz,  [4].  The 
first  real  data  represent  three  corner  reflectors  rotating  at  approximately  60  RPM 
(rotation  per  minute)  and  the  rigid  body  observed  by  the  radar  with  Tr  =  1  kHz. 
The  STFT  of  the  returned  signal,  for  the  given  range  bin,  is  shown  in  Figure  20(a). 
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Figure  20:  Real  radar  data  corresponding  to  a  rigid  body  and  three  corner  reflectors 
rotating  at  ~  60  RPM  (a -d).  (a)  the  STFT.  (b)  Sorted  STFT.  (c)  The  original  FT 
and  (d)  the  FT  reconstructed  by  summing  over  50%  of  the  lowest  STFT  samples. 
The  same  procedure  is  repeated  for  real  radar  data  corresponding  to  a  stronger  rigid 
body  and  two  corner  reflectors  rotating  at  ~  40  RPM  (e-h).  A  logarithmic  amplitude 
scale  is  used  in  subplots  (g)  and  (h). 
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After  sorting  the  STFT  over  time  Figure  20(b),  the  constant  frequency  component 
corresponding  to  the  rigid  body  becomes  more  visible,  since  the  time  varying  fre¬ 
quency  content  is  spread  over  many  frequencies,  for  each  frequency  bin.  The  rigid 
body  is  separated  from  the  m-D  and  its  FT  is  successfully  reconstructed  by  using 
50%  of  the  lowest  STFT  values,  as  shown  in  Figure  20(d).  If  we  compare  it  to  the 
FT  of  the  original  signal,  Figure  20(c),  we  can  see  the  improvement  in  the  rigid  body 
presentation. 

In  the  second  example,  the  real  radar  data  corresponding  to  two  outside  corner  re¬ 
flectors,  rotating  at  approximately  40  RPM  (all  facing  radar)  with  rigid  body,  are 
analyzed.  The  same  radar  as  in  the  previous  example  is  used,  while  the  reflectivity 
of  rigid  body  is  much  higher  than  those  of  the  rotating  reflectors.  The  STFT  repres¬ 
entation  of  the  observed  signal  is  shown  in  Figure  20(e).  The  sorted  STFT  is  shown 
in  Figure  20(f).  The  original  FT  is  shown  in  Figure  20(g).  The  reconstructed  FT, 
obtained  by  summing  50%  of  the  lowest  STFT  values  is  presented  in  Figure  20(h). 
We  have  successfully  removed  most  of  the  m-D.  Moreover,  we  may  use  the  removed 
STFT  samples  in  order  to  estimate  features  of  rotating  reflectors.  Here,  we  have  used 
a  logarithmic  scale  to  present  the  reconstructed  values,  since  the  m-D  values  were 
very  low. 

4  Micro-Doppler  toolbox 


The  micro-Doppler  toolbox  is  a  set  of  functions  for  the  analysis  of  m-D  in  radar  return 
signals.  Two  approaches  are  implemented. 

•  L-statistics  is  used  in  order  to  separate  the  m-D  and  the  rigid  body  from  the 
analysed  range  bin. 

•  Micro-Doppler  parameter  estimation  is  performed  in  order  to  estimate  the  m-D 
frequency,  amplitude  and  phase.  The  parameter  estimation  method  based  on 
the  inverse  Radon  transform  is  developed.  The  developed  method  is  capable  to 
detect  m-D  parameters  even  in  the  case  of  multicomponent  m-D  signal. 

•  The  toolbox  is  capable  to  deal  with  one-dimensional  (ID)  and  two-dimensional 
(2D)  signals.  If  the  2D  signal  is  analyzed,  a  procedure  for  detection  range  bins 
with  possible  m-D  is  implemented. 

•  The  toolbox  offers  link  to  the  Time-Frequency  Analysis  toolbox  in  order  to 
provide  highly  concentrated  time-frequency  representation  of  the  analyzed  sig¬ 
nal. 

Beside  core  functions  for  L-statistics  and  inverse  Radon  transform-based  parameter 
estimation,  the  toolbox  contains  Virtual  Instrument  and  Demo. 
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4.1  Files  and  installation 

This  toolbox  has  six  program  hies: 


L  statist  ics.m  The  function  takes  single  range  bin  signal  (in  time  domain)  and 
performs  separation  of  rigid-body  and  m-D  parts  of  the  analyzed  signal.  For 
more  information  type  help  L_statistics  in  command  line. 

MD  iradon.m  This  function  estimates  parameters  of  m-D  component.  It  is  suit¬ 
able  for  single-component  m-D  signals.  In  multi-component  signals,  the  m-D 
function  estimates  parameters  for  one  component  only. 

MD  iradonMC . m  This  function  estimates  parameters  of  m-D  in  multi-component 
case.  After  the  estimation  of  parameters,  the  estimated  component  is  removed 
from  the  analyzed  signal  and  the  estimation  procedure  is  repeated  for  the  next 
component. 

MD  VI.rn  Virtual  Instrument  for  m-D  analysis. 

MD  demo.m  Demonstration  of  core  functions  usage. 

STFT  TFR.m  Short-time  Fourier  transform  calculation.  This  function  is  part  of 
the  Time-Frequency  Analysis  toolbox. 

SM_TFR.m  S-method  calculation.  This  function  is  part  of  the  Time-frequency 
Analysis  toolbox. 

test_mat  file. mat  Matlab  hie  with  ID  and  2D  data  variables.  This  hie  is  used 
for  testing  the  External  data  ID  and  External  data  2D  functionality. 

To  install  the  toolbox,  copy  all  hies  in  the  working  directory.  Alternatively,  hies  can 

be  copied  in  the  directory  that  includes  MATLAB  path. 

4.2  Virtual  instrument 

•  Virtual  Instrument  is  started  with  command  MD_VI.  The  main  screen  is  dis¬ 
played  in  Figure  21. 


4.2.1  Main  screen 

•  The  main  screen  is  divided  into  three  parts.  The  left  part  is  dedicated  to  the 
definition  of  the  input  data.  The  center  part  deals  with  L-statistics  as  a  method 
for  separating  m-D  and  rigid  body  parts  of  the  analyzed  signal.  The  right  part 
is  for  the  m-D  parameter  estimation. 
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Figure  21:  Virtual  Instrument  -  main  screen 
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Figure  22:  Load  ID  external  data 
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•  The  input  data  can  be  simulated  signals,  ID  external  data  or  2D  external  data. 

•  If  example  signals  are  selected,  the  user  can  choose  one  of  six  predefined  data 
sets,  where  Example  signals  1-4  and  Example  signal  6  are  ID  data  and  Example 
signal  5  is  2D  data. 

•  When  the  example  signal  is  selected,  the  ToolTip  with  explanation  of  the  chosen 
signal  is  displayed  if  the  user  hold  the  mouse  over  the  example  list. 

•  In  the  2D  data  case,  the  edit  box  with  range-bin  selection  and  the  button  for 
search  of  m-D  range  bins  are  enabled. 

•  When  the  user  select  the  data,  it  is  allowed  to  plot  TF  representation  of  the 
selected  data. 

•  The  pulse  repetition  frequency  can  be  defined  in  the  lower  part  of  the  left  side 
of  the  Virtual  Instrument. 

•  If  ID  external  data  is  selected,  the  user  can  select  the  data  format  (time  do¬ 
main  or  frequency  domain)  and  load  the  data  from  external  MATLAB  hie  as 
presented  in  Figure  22. 

•  If  variable  q  is  detected  in  the  input  hie,  q  is  used  as  input  data.  Otherwise, 
the  user  can  choose  ID  variable  from  the  data  hie  to  use  it  as  input  data. 

•  When  2D  external  data  is  selected  (Figure  23),  the  user  has  choices  as  in  the 
previous  case  except  the  user  can  now  transpose  the  input  data  if  the  range  bins 
are  not  in  rows,  and  the  range  bin  should  be  selected  prior  to  further  analysis. 

•  In  this  case,  if  the  variable  q  is  not  found  in  the  selected  MAT  hie,  the  user  can 
select  2D  variable  that  contains  radar  data. 
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Figure  23:  Load  2D  external  data 

4.2.2  Time-frequency  data  representation 

•  Click  on  the  TF  representation  button  that  gives  the  TF  representation  of  item 
the  selected  data.  Example  4  is  presented  in  Figure  24.  The  S-method  with 
L  =  2  and  the  appropriate  time  domain  window  is  used  in  this  example. 

4.2.3  L-statistics 

•  The  middle  part  of  the  Virtual  Instrument  screen  (Figure  21)  is  dedicated  to 
L-statistics. 

•  The  user  can  select  the  window  type  and  window  length  (or  determine  the 
window  length  automatically  according  to  concentration  measures). 

•  In  the  second  step,  the  user  can  define  thresholds  for  the  m-D  and  rigid  body 
parts  separation  or  can  choose  adaptive  threshold. 

•  In  the  third  step,  the  user  can  define  the  time-step  for  the  STFT  or  use  auto¬ 
matically  determined  time-step. 
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•  Click  on  the  Analyze  data  button  provides  results  in  a  form  similar  to  Figure 
25. 

•  The  upper  left  subplot  is  the  TFR  (time-frequency  representation)  of  the  ana¬ 
lyzed  data.  The  sorted  values  of  the  TFR  along  with  threshold(s)  are  presented 
in  the  upper  right  subplot. 

•  The  Fourier  transform  of  the  original  signal  and  the  rigid  body  part  are  presen¬ 
ted  in  the  lower  left  two  subplots. 

•  The  lower  right  subplot  is  the  TFR  of  the  m-D  part  of  the  analyzed  signal. 


Figure  24:  Time- frequency  representation  of  the  analyzed  data 

4.2.4  Micro-Doppler  parameter  estimation 

•  The  third  part  of  the  Virtual  Instrument  is  dedicated  to  the  m-D  parameter 
estimation  based  on  the  inverse  Radon  transform. 

•  The  user  can  select  between  the  single  or  multi  component  case. 

•  Furthermore,  the  user  can  give  expected  number  of  periods  within  the  analyzed 
image  or  define  the  period  range. 

•  In  cases  when  the  m-D  signal  is  not  centered  around  the  middle  cross-range  bin, 
the  user  can  define  the  cross-range  shift  or  chose  the  shift  automatically. 
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•  Click  on  the  Estimate  m-D  button  when  single  component  m-D  is  selected;  this 
gives  results  as  illustrated  in  Figure  26  and  Figure  27. 

•  Figure  25  presents  the  TF  representation  of  the  analyzed  data  along  with  es¬ 
timated  m-D  (black  line). 

•  Estimated  m-D  parameters  (Rate,  Amplitude,  Phase  and  the  number  of  peri¬ 
ods)  are  given  in  the  message-box  presented  in  Figure  27. 


Figure  25:  L-statistics  analysis 
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Figure  26:  Micro-Doppler  parameter  estimation  results 


)  Estimated  pa...  -  □ 


E  stimated  parameters: 


m-D  rate  [1/s] :  0.97217 
m-D  normalized  amplitude  :  0.85819 
m-D  initial  phase  [deg] :  -119.5785 
Number  of  periods  :  1.001 6 


OK 


Figure  27:  Micro-Doppler  parameters 
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Figure  28:  Micro-Doppler  parameter  estimation  results  -  multicomponent  case 
second  analyzed  component 
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Figure  29:  Micro-Doppler  parameters  for  second  component  in  multi-component  case 


When  the  multicomponent  m-D  is  selected,  the  search  procedure  is  repeated  in  order 
to  detect  all  m-D  components.  Results  are  presented  in  Figure  28  for  the  second  com¬ 
ponent  in  the  analyzed  signal.  Parameters  of  the  analyzed  component  are  displayed 
and  the  user  can  choose  to  search  for  another  component  or  to  stop  analysis  (Figure 
29). 

4.2.5  2D  data  analysis 

•  If  2D  data  is  analyzed,  the  user  can  activate  the  search  for  m-D  bins  resulting 
in  Figure  30.  Here  the  user  can  select  the  detected  range  bin  from  the  list  on 
the  right  and  go  back  to  the  Virtual  Instrument  to  analyze  the  selected  range 
bin. 

4.2.6  Demo 

•  The  usage  of  core  functions  L_statistics,  MD  iradon  and  MD  irdonMC  is 
demonstrated  with  demo  program  as  presented  in  Figure  31. 

4.3  Time-Frequency  Toolbox 

The  time-frequency  analysis  toolbox  is  a  set  of  basic  functions  used  in  other  toolboxes 
such  as  SAR/ISAR  toolbox  and  micro-Doppler  toolbox.  However,  the  usage  of  this 
toolbox  is  not  limited  to  the  support  to  other  toolboxes.  One  can  use  this  toolbox  in 
order  to  get  the  time-frequency  representation  of  an  arbitrary  input  signal. 
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Figure  30:  Range  bin  selection  with  2D  input  data 


Figure  31:  Range  bin  selection  with  2D  input  data 
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The  toolbox  provides  two  time- frequency  methods:  short-time  Fourier  transform 
(STFT)  and  S-method  (SM).  A  variety  of  options  are  added  to  the  core  functions, 
allowing  the  user,  for  example,  to  select  the  window  with  constant  width  or  auto¬ 
matically  determine  the  window  width  in  STFT  method  and  use  constant  $L$  or 
adaptive  S-method.  Most  of  the  options  are  demonstrated  by  example  signals  with 
hie  TFR_demo .  m. 

4.3.1  Files  and  installation 

The  toolbox  has  four  program  hies: 

STFT_TFR.m 
SM  TFR.m 
ShowData.m 
TFRdemo.m 

To  install  the  toolbox,  copy  all  the  hies  in  the  working  directory. 


4.3.2  Demo 

The  main  toolbox  window  is  presented  in  Figure  32.  The  user  can  select  an  example 
signal  by  using  the  list  on  the  left.  Thereafter  the  appropriate  STFT  command  should 
be  chosen  and  finally  the  S-method  command  should  be  selected  on  the  right  part. 

The  user  can  now  click  for  STFT  or  S-method  calculation  and  obtain  the  results  as 
presented  in  Figure  33  for  STFT  and  Figure  34  for  S-method. 


4.3.3  Advanced  visualization 

The  advanced  visualization  options  are  invoked  by  typing  ShowData(SM)  at  the  MAT- 
LAB  command  line  for  S-method  or  ShowData(ST)  for  STFT.  The  result  is  presented 
in  Figure  35. 
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Figure  32:  Time- frequency  analysis  toolbox  demo. 
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Figure  33:  Short-Time  Fourier  Transform  of  the  analyzed  signal. 
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Figure  34:  S-method  of  the  analyzed  signal. 
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Figure  35:  Advanced  visualization  options. 
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5  Conclusion 


A  method  for  estimation  of  the  parameters  of  sinusoidally  modulated  signal  is  in¬ 
troduced.  The  proposed  method  is  based  on  the  inverse  Radon  transform  and  the 
concentration  measures.  It  is  shown  that  proposed  method  provides  promising  estim¬ 
ation  and  decomposition  results  for  monocomponent  and  multicomponent  signals. 
The  noise  and  interferences  influence  upon  the  estimation  procedure  is  considered. 
It  can  be  concluded  that  the  proposed  method  is  very  robust  to  the  noise  and  other 
interferences.  We  have  also  shown  that  the  results  obtained  by  the  proposed  method 
are  meaningful  even  in  cases  when  the  analyzed  signal  is  periodic  but  not  sinusoid¬ 
ally  modulated.  It  can  be  used  to  estimate  the  parameters  of  periodic  extension  of  a 
non-periodic  time-frequency  patterns  and  partially  available  data  as  well. 

The  m-D  effect  appears  in  the  ISAR/SAR  images  in  the  case  when  there  are  fast 
moving  reflectors  in  the  observed  scene.  The  m-D  can  severely  decrease  the  quality 
and  readability  of  the  obtained  radar  image.  Its  detection  and  removal  are  very 
important  for  obtaining  a  focused  image  of  the  rigid  body.  In  the  report,  an  algorithm 
for  the  m-D  removal  and  reconstruction  of  the  rigid  body  image  are  proposed.  It  is 
based  on  the  L-statistics.  Although  the  L-statistics  are  already  applied  by  the  authors 
for  a  similar  purpose,  the  algorithm  proposed  here  is  simpler  and  produces  better 
results.  The  reconstruction  is  performed  by  using  only  one  complex  STFT,  rather 
than  using  several  STFTs  (with  absolute  values)  and  order  statistics  combination, 
as  it  was  done  in  the  previous  work.  Since  the  proposed  algorithm  is  based  on 
the  L-statistics,  being  a  tool  for  robust  signal  analysis,  it  is  robust  to  the  effects  of 
noise.  Also,  it  can  behave  as  a  high-resolution  algorithm,  since  it  can  separate  close 
rigid  body  points.  In  order  to  improve  the  computational  efficiency,  an  adaptive 
threshold  is  used  to  distinguish  among  STFT  samples  that  correspond  to  the  rigid 
body  and  moving  parts.  Moreover,  two  additional  thresholds  are  incorporated  in  the 
algorithm.  The  first  threshold  detects  whether  there  is  a  returned  radar  signal  in  a 
range  bin,  while  the  second  threshold  detects  whether  there  exist  m-D  effects  in  a 
range  bin.  Consequently,  the  procedure  for  the  m-D  extraction  could  be  performed 
only  for  the  range  bins  where  the  m-D  effect  is  detected;  that  could  lead  to  the 
overall  computational  savings.  Through  the  examples,  it  is  shown  that  the  proposed 
algorithm  successfully  separates  the  rigid  body  and  the  m-D  effects. 

Based  on  L-statistics,  spectrogram,  and  inverse  Radon  transform,  a  radar  data  ex¬ 
ploitation  Matlab  toolbox  was  developed  for  classifying  air  and  land  targets.  This 
toolbox  will  scan  and  search  any  M-D  activity  in  the  desired  area  of  the  SAR/ISAR 
image.  The  toolbox  will  have  three  options;  1)  Focus  the  stationary  target  if  there 
is  no  M-D,  2)  Focus  the  target  after  removing  the  M-D,  and  3)  extract  the  M-D 
and  determine  the  motion  parameters  such  as  rotation/ vibration  rate,  initial  phase, 
and  amplitude  of  the  stationary  and  moving  targets.  The  toolbox  was  tested  and 
validated  against  a  variety  of  simulated  and  measured  targets. 


60 


DRDC-RDDC-201 4-R6 


References 


[1]  Thayaparan,  T.,  Abrol,  S.,  Riseborough,  E.,  Stankovic,  L.,  Lamothe,  D.,  and 
Duff,  G.  (2007).  Analysis  of  radar  micro- Doppler  signatures  from  experimental 
helicopter  and  human  data,  IET  Proceedings  Radar  Sonar  Navig .,  Vol.  1,  No.  4, 
pp.  288-299. 

[2]  Thayaparan,  T.,  Suresh,  P.,  and  Qian,  S.  (2010).  Micro-Doppler  Analysis  of 
Rotating  Target  in  SAR,  IET  Signal  Processing,  IET  Signal  Processing ,  Vol.  4, 
No.  3,  pp.  245-255. 

[3]  Thayaparan,  T.,  Stankovic,  L.,  Dakovic,  M.,  and  Popovic,  V.  (2010), 
Micro-Doppler  parameter  estimation  from  a  fraction  of  the  period,  IET  Signal 
Processing ,  Vol.  4,  No.  3,  pp.  201-212. 

[4]  Thayaparan,  T.,  Stankovic,  L.,  and  Djurovic,  I.  (2008).  Micro-Doppler  Based 
Target  Detection  and  Feature  Extraction  in  Indoor  and  Outdoor 
Environments,  J.  of  the  Franklin  Institute ,  Vol.  345,  No.  6,  pp.  700-722. 

[5]  Stankovic,  L.,  Thayaparan,  T.,  and  Djurovic,  I.  (2006).  Separation  of  target 
rigid  body  and  micro-Doppler  effects  in  ISAR  imaging,  IEEE  Trans.  Aerosp. 
Electron.  Syst.,  Vol.  41,  No.  4.  pp.  1,496-1,506. 

[6]  Chen,  V.  C.,  Li,  F.,  Ho,  S.-S.,  and  Wechsler,  H.  (2006).  Analysis  of 
micro-Doppler  signatures,  IEE  Proc.  Radar,  Sonar,  Navig.,  Vol.  150,  No.  4,  pp. 
271-276. 

[7]  Chen,  V.  C.  (2006).  Micro-Doppler  effect  in  radar:  Part  I:  Phenomenon, 
physics,  mathematics,  and  simulation  study,  IEEE  Trans,  on  Aerosp.  Electron. 
Syst ,  Vol.  42,  No.  1. 

[8]  Stankovic,  L.,  Dakovic,  M,  and  Thayaparan,  T.  (2013).  Time- Frequency  Signal 
Analysis  With  Applications,  Artech  House. 

[9]  Chen,  V.  C.  and  Ling,  H.  (2002).  Time- frequency  transforms  for  radar  imaging 
and  signal  analysis,  Artech  House,  Boston. 

[10]  Li,  J.  and  Ling,  H.  (2003).  Application  of  adaptive  chirplet  representation  for 
ISAR  feature  extraction  from  targets  with  rotating  parts,  IEE  Proc.  Radar, 
Sonar,  Navig.,  Vol.  150,  No.  4,  pp.  284-291. 


DRDC-RDDC-201 4-R6 


61 


[11]  Sabry,  R.,  Vachon,  P.  W.,  Sandirasegaram,  N.,  Seeker,  J.,  Liu,  C.,  Mattar,  K. 
E.,  Wong,  S.,  and  Thayaparan,  T.  (2013).  Space-based  Radar  Products  for 
Military  Operations  (U):  Project  15EL  Final  Report  (U),  DRDC  Ottawa 
Research  Centre,  TR-2013-125. 

[12]  Mattar,  K.E.,  Guertin,  R.,  Schlingmeier,  D.,  Thayaparan,  T.,  Lukowski,  T.I, 
Sabry,  R.,  Seeker,  J.,  Barrie,  G.  (2013).  Makeshift  radar  reflector  trial  for 
search  and  rescue  scenarios.  DRDC  Ottawa  Research  Centre,  TM-2013-084. 

[13]  Mattar,  K.E.,  Barrie,  G.,  English,  R.A.,  Liu,  C.,  Sabry,  R.,  Sandirasegaram, 

N.,  Seeker,  J.,  Thayaparan,  T.  And,  Wolfe,  J.  (2011).  An  assessment  of 
Synthetic  Aperture  Radar  feature  extraction  applications,  Ottawa  TM 
2011-227,  DRDC  Ottawa. 

[14]  Vachon,  P.W.,  Barrie,  G.,  Baziuk,  W.,  English,  R.,  Liu,  C.,  Mattar,  K.,  Sabry, 
R.,  Sandirasegaram,  N.,  Schlingmeier,  D.,  Seeker,  J.,  Thayaparan,  T.,  Wilcox., 
C.,  (2009).  Synthetic  Aperture  Radar  exploitation:  RDE  Groups’  perspective, 
DRDC  Ottawa,  TR  2009-194. 

[15]  Thayaparan,  T.,  Stankovic,  LJ,  and  Djorovic,  I.  (2007).  Target  Detection  and 
Feature  Extraction  in  Indoor  and  Outdoor  Environments  using  micro-Doppler 
analysis,  DRDC  Ottawa,  TM  2008-255. 

[16]  Thayaparan,  T.  (2006)  Separation  of  target  rigid  body  and  micro-Doppler 
effects  in  ISAR/SAR,  DRDC  Ottawa,  TM  2006-187. 

[17]  Thayaparan,  T.  and  Abrol,  S.  (2005)  Micro-Doppler  Analysis  of  Rotating 
Target  in  SAR,  DRDC  Ottawa  TM  2005-204,  December  2005. 

[18]  Thayaparan,  T.,  Abrol,  S.  and  Riseborough,  E.  (2004).  Micro-Doppler  radar 
signatures  for  intelligent  target  recognition,  DRDC  Ottawa  TM  2004-170. 

[19]  P.  Suresh,  T.  Thayaparan,  K.  Venkataramaniah,  Micro-Doppler  Analysis  of 
Rotating  Targets  using  the  Adaptive  S-  method,  International  Radar 
Symposium  India  2011,  30  Nov-4  Dec  2011,  Bangalore,  India. 

[20]  Thayaparan,  T.,  Stankovic,  L.,  Darkovic,  M.,  Suresh,  P.,  Venkataramaniah,  K., 
SivaSankaraSai,  S.,  Shankar,  S.,  Sairam,  T.,  Nikhilesh.  K.  (2010).  Parameter 
estimation  from  a  fraction  of  the  period  for  a  rotating  target,  International 
Radar  Symposium  2010,  June  16-18,  Lithuania. 

[21]  Suresh,  P.,  Thayaparan,  T.,  Venkataramaniah,  K.,  SivaSankaraSai,  S.,  and 
Sridharan,  K.  S.  (2009).  Analysis  of  micro- Doppler  radar  signatures  in  SAR 
using  S-method-based,  International  Radar  Symposium  India-  2009,  December 
8-11,  Bangalore,  India. 


62 


DRDC-RDDC-201 4-R6 


[22]  Suresh,  P.,  Thayaparan,  T.,  Venkataramaniah,  K.,  SivaSankaraSai,  S.,  and 
Sridharan,  K.  S.  (2009).  Extracting  micro  Doppler  radar  signatures  from 
experimental  helicopter  data  using  adaptive  Chirplet-based  analysis, 
International  Radar  Symposium  India-  2009,  December  8-11,  Bangalore,  India. 

[23]  Thayaparan,  T.,  Dakovic,  M.,  Stankovic,  L.,  and  Suresh,  P.  (2009).  Intelligent 
target  recognition  using  micro-Doppler  radar  signatures,  2009  SPIE  Defense, 
Security,  and  Sensing,  Proceedings  of  SPIE  Volume  7308. 

[24]  Thayaparan,  T.,  Stankovic,  L.,  and  Dakovic,  M.  (2007).  Micro-Doppler  Feature 
Extraction  Using  S-method,  International  Radar  Symposium  -  IRS  2006,  5-7 
September,  Colonge,  Germany,  pp.  833-837. 

[25]  Thayaparan,  T.,  Micro-Doppler  Feature  Extraction  using  Time-Frequency 
Analysis,  NATO  SET-155,  Advancing  Sensing  Thorough  the  Wall  (STTW) 
Technologies,  October  5-7,  2009 

[26]  Thayaparan,  T.  ,  Sensing  and  Communications  Using  Ultra-wideband  Random 
Noise  Radar,  TTCP  MAR  AG-10  -  Maritime  Force  Protection,  CFMWC, 
MARLANT,  Halifax,  Nova  Scotia,  Canada  April  20-24,  2009. 

[27]  Thayaparan,  T.  and  Abrol,  Chen,  V.  C.,  Analysis  of  micro-Doppler  radar 
signatures  from  experimental  helicopter  and  human  data,  NATO  SET-080 
symposium,  Oslo,  Norway,  October  11-13,  2004. 

[28]  Thayaparan,  T.,  Non-stationary  Radar  Signal  Analysis:  Time-Frequency 
Approach,  delivered  3-hour  Tutorial,  International  Radar  Symposium  India 
2011,  30  Nov-4  Dec  2011,  Bangalore,  India 

[29]  Barbarossa,  S.  (1995).  Analysis  of  multicomponent  LFM  signals  by  a  combined 
Wigner-Hough  transform,  IEEE  Trans,  on  Signal  Processing,  Vol.  43,  No.  6, 

pp.  1,511-1,515. 

[30]  Barbarosa  S.  and  Lemoine,  O.  (1996).  Analysis  of  nonlinear  FM  signals  by 
pattern  recognition  of  their  time-frequency  representation,  IEEE  Signal 
Processing  Letters,  Vol.  3,  No.  4,  pp.  112-115. 

[31]  Stankovic,  L.,  Thayaparan,  T.,  Dakovic,  M.,  and  Popovic-Bugarin,  V.  (2013). 
Micro-Doppler  removal  in  the  radar  imaging  analysis,  IEEE  Trans,  on 
Aerospace  and  Electronics  Systems,  Vol.  49,  No.  2,  pp.  1,234-1,250. 

[32]  Stankovic,  L.,  Thayaparan,  T.,  and  Djurovic,  I.  (2006).  Separation  of  Target 
Rigid  Body  and  Micro-Doppler  Effects  in  ISAR  Imaging,  IEEE  Trans,  on 
Aerospace  and  Electronics  Systems,  Vol.  42,  No.  7,  pp.  1,496-1,506. 


DRDC-RDDC-201 4-R6 


63 


[33]  Ristic,  B.  and  Boashash,  B.  (1993).  Kernel  design  for  time-frequency  signal 
analysis  using  the  Radon  transform,  IEEE  Trans,  on  Signal  Processing,  Vol.  41, 
No.  5,  pp.  1,996-2,008. 

[34]  Stankovic,  S.,  Djurovic,  I.,  and  Pitas,  I.  (2001).  Watermarking  in  the  space 
/spatial- frequency  domain  using  two-dimensional  Radon-Wigner  distribution, 
IEEE  Trans,  on  Image  Processing ,  Vol.  10,  No.  4,  pp.  650-658. 

[35]  Stankovic,  L.  (2001).  A  measure  of  some  time-frequency  distributions 
concentration,  Signal  Processing,  Vol.  81,  No.  3,  pp.  621-631. 

[36]  Wood,  J.  C.  and  Barry,  D.  T.  (1994).  Linear  signal  synthesis  using  the 
Radon-Wigner  distribution,  IEEE  Trans,  on  Signal  Processing,  Vol.  42,  No.  8, 
pp.  2,105-2,111. 

[37]  Wood,  J.  C.  and  Barry,  D.  T.  (1994).  Radon  transformation  of  time- frequency 
distributions  for  analysis  of  multicomponent  signals,  IEEE  Trans,  on  Signal 
Processing,  Vol.  42,  No.  11,  pp.  3,166-3,177. 

[38]  Herman,  G.  T.  (1980).  Image  Reconstruction  from  Projections:  The 
Fundamentals  of  Computerized  Tomography ,  Academic  Press,  New  York. 

[39]  Sommer,  H.  and  Salerno,  J.  (1971).  Radar  target  identification  system,  U.S. 
Patent  3,  pp.  614-779. 

[40]  Marple,  S.  L.  (2001).  Large  dynamics  range  time- frequency  signal  analysis  with 
application  to  helicopter  Doppler  radar  data,  Sixth  International  Symposium 
on  Signal  Processing  and  its  Applications  -  ISSPA,  13-16  August,  Kuala 
Lumpur,  Malaysia. 

[41]  Zediker,  M.  S.,  Rice,  R.  R.,  and  Hollister,  J.  H.  (1998).  Method  for  extending 
range  and  sensitivity  of  a  fiber  optic  nricro-Doppler  system  and  apparatus,  U.S. 
Patent  6,847,817. 

[42]  Stankovic,  L.,  Dakovic,  M.,  and  Ivanovic,  V.  N.  (2001).  Performance  of 
spectrogram  as  IF  estimator,  Electronics  Letters,  Vol.  37,  No.  12,  pp.  797-799. 

[43]  Stankovic,  L.  (2004).  Performance  Analysis  of  the  Adaptive  Algorithm  for 
Bias-to- Variance  Trade-off,  IEEE  Trans,  on  Signal  Processing,  Vol.  52,  No.  5, 
pp.  1,228-1,234. 

[44]  Bai,  X.,  Zhou,  F.,  Xing,  M.,  and  Bao,  Z.  (2011).  High  resolution  ISAR  imaging 
of  targets  with  rotating  parts,  IEEE  Trans,  on  Aerosp.  Electron.  Syst,  Vol.  47, 
No.  4,  pp.  2,530-2,543. 

[45]  Totir,  F.  and  Radoi,  E.  (2009).  Superresolution  algorithms  for  spatial  extended 
scattering  centers,  Digital  Signal  Processing,  Vol.  19,  No.  5,  pp.  780-792. 


64 


DRDC-RDDC-201 4-R6 


[46]  Martorella,  M.  (2008).  Novel  approach  for  ISAR  image  cross-range  scaling, 
IEEE  Trans.  Aerosp.  Electron.  Syst .,  Vol.  44,  No.  1,  pp.  281-294. 

[47]  Martorella,  M.  and  Berizzi,  F.  (2005).  Time  windowing  for  highly  focused 
ISAR  image  reconstruction,  IEEE  Trans.  Aerosp.  Electron.  Syst.,  Vol.  41,  No. 
3,  pp.  992-1,007. 

[48]  Wang,  Y.  and  Jiang,  Y.-C.  (2011).  ISAR  imaging  of  ship  target  with  complex 
motion  based  on  new  approach  of  parameters  estimation  for  polynomial  phase 
signal,  EURASIP  Journal  on  Advances  in  Signal  Processing ,  Volume  2011, 
Article  ID  425203. 

[49]  Sparr,  T.  and  Krane,  B.  (2003).  Micro-Doppler  analysis  of  vibrating  targets  in 
SAR,  IEE  Proc.  Radar  Sonar  Navig.,  Vol.  150,  No.  4,  pp.  277-283. 

[50]  Lyonnet,  B.,  Ioana,  C.,  and  Amin,  M.  G.  (2010).  Human  gait  classification 
using  microDoppler  time-frequency  signal  representations,  Radar  Conference 
2010  IEEE ,  10-14  May,  Washington,  DC,  pp.  915-919. 

[51]  Djurovic,  I.,  Stankovic,  L.,  Bohme,  J.  F.  (2003).  Robust  L-estimation  based 
forms  of  signal  transforms  and  time-frequency  representations,  IEEE  Trans. 
Signal  Processing ,  Vol.  51,  No.  7,  pp.  1,753  -  1,761. 

[52]  Wang,  Y.,  Ling,  H.,  and  Chen,  V.  C.  (1998).  ISAR  motion  compensation  via 
adaptive  joint  time- frequency  techniques,  IEEE  Trans.  Aerosp.  Electron.  Syst., 
Vol.  38,  No.  2,  pp.  670-677. 

[53]  Katkovnik,  V.  and  Stankovic,  L.  (1998).  Instantaneous  frequency  estimation 
using  the  Wigner  distribution  with  varying  and  data  driven  window  length, 
IEEE  Trans.  Signal  Processing ,  Vol.  46,  No.  9,  pp.  2,315-2,325. 

[54]  Yin,  P.-Y.  (1999).  A  new  circle/ellipse  detector  using  genetic  algorithms, 
Pattern  Recognition  Letters,  Vol.  20,  No.  7,  pp.  731-740. 

[55]  Misiurewicz,  J.,  Kulpa,  K.,  and  Czekala,  Z.  (1998).  Analysis  of  recorded 
helicopter  echo,  Radar  97,  14-16  Oct,  Edinburgh,  UK.,  pp.  449-453. 

[56]  Marple,  S.  L.,  (2004).  Special  time-frequency  analysis  of  helicopter  Doppler 
radar  data,  in  Time- Frequency  Signal  Analysis  and  Processing ,  ed.  B. 
Boashash,  Elsevier. 

[57]  Djurovic,  I.,  Thayaparan,  T.,  and  Stankovic,  L.  (2006).  Adaptive  local 
polynomial  Fourier  transform  in  ISAR,  EURASIP  Journal  of  Applied  Signal 
Processing,  Vol.  2006,  Article  ID  36093. 

[58]  Stankovic,  L.  (2001).  A  measure  of  some  time-frequency  distributions 
concentration,  Signal  Processing,  Vol.  81,  No.  3,  pp.  621-631. 


DRDC-RDDC-201 4-R6 


65 


This  page  intentionall  left  blank. 


66 


DRDC-RDDC-201 4-R6 


DOCUMENT  CONTROL  DATA 

(Security  markings  for  the  title,  abstract  and  indexing  annotation  must  be  entered  when  the  document  is  Classified  or  Designated.) 

2a.  SECURITY  MARKING  (Overall  security 
marking  of  the  document,  including 
supplemental  markings  if  applicable.) 

UNCLASSIFIED 

2b.  CONTROLLED  GOODS 

(NON-CONTROLLED 
GOODS) 

DMCA 

REVIEW:  GCEC  APRIL  2011 

3.  TITLE  (The  complete  document  title  as  indicated  on  the  title  page.  Its  classification  should  be  indicated  by  the  appropriate 
abbreviation  (S,  C  or  U)  in  parentheses  after  the  title.) 

A  Virtual  Instrument  for  Micro-Doppler  Analysis  of  Signals  in  SAR/ISAR:  Parameters 
Estimation  and  Focusing  Moving  Targets 

4.  AUTHORS  (Last  name,  followed  by  initials  -  ranks,  titles,  etc.  not  to  be  used.) 

Thayaparan,  T. 


1 .  ORIGINATOR  (The  name  and  address  of  the  organization  preparing  the 
document.  Organizations  for  whom  the  document  was  prepared,  e.g.  Centre 
sponsoring  a  contractor's  report,  or  tasking  agency,  are  entered  in  section  8.) 

DRDC  -  Ottawa  Research  Centre 

3701  Carling  Avenue,  Ottawa  ON  K1 A  0Z4,  Canada 


5.  DATE  OF  PUBLICATION  (Month  and  year  of  publication  of 
document.) 

6a.  NO.  OF  PAGES  (Total 
containing  information. 

Include  Annexes, 

Appendices,  etc.) 

6b.  NO.  OF  REFS  (Total 
cited  in  document.) 

April  2014 

80 

58 

7.  DESCRIPTIVE  NOTES  (The  category  of  the  document,  e.g.  technical  report,  technical  note  or  memorandum.  If  appropriate,  enter 
the  type  of  report,  e.g.  interim,  progress,  summary,  annual  or  final.  Give  the  inclusive  dates  when  a  specific  reporting  period  is 
covered.) 


Scientific  Report 


8.  SPONSORING  ACTIVITY  (The  name  of  the  department  project  office  or  laboratory  sponsoring  the  research  and  development  - 
include  address.) 

DRDC  -  Ottawa  Research  Centre 

3701  Carling  Avenue,  Ottawa  ON  K1 A  0Z4,  Canada 


9a.  PROJECT  OR  GRANT  NO.  (If  appropriate,  the  applicable 
research  and  development  project  or  grant  number  under 
which  the  document  was  written.  Please  specify  whether 
project  or  grant.) 

9b.  CONTRACT  NO.  (If  appropriate,  the  applicable  number  under 
which  the  document  was  written.) 

15dp 

f  0a.  ORIGINATOR’S  DOCUMENT  NUMBER  (The  official 

document  number  by  which  the  document  is  identified  by  the 
originating  activity.  This  number  must  be  unique  to  this 
document.) 

f  0b.  OTHER  DOCUMENT  NO(s).  (Any  other  numbers  which  may 
be  assigned  this  document  either  by  the  originator  or  by  the 
sponsor.) 

DRDC-RDDC-201 4-R6 

1 1 .  DOCUMENT  AVAILABILITY  (Any  limitations  on  further  dissemination  of  the  document,  other  than  those  imposed  by  security 
classification.) 


(X)  Unlimited  distribution 

(  )  Defence  departments  and  defence  contractors;  further  distribution  only  as  approved 
(  )  Defence  departments  and  Canadian  defence  contractors;  further  distribution  only  as  approved 
(  )  Government  departments  and  agencies;  further  distribution  only  as  approved 
(  )  Defence  departments;  further  distribution  only  as  approved 
(  )  Other  (please  specify): 


12.  DOCUMENT  ANNOUNCEMENT  (Any  limitation  to  the  bibliographic  announcement  of  this  document.  This  will  normally  correspond 
to  the  Document  Availability  (11).  However,  where  further  distribution  (beyond  the  audience  specified  in  (11))  is  possible,  a  wider 
announcement  audience  may  be  selected.)  UNLIMITED 


1 3.  ABSTRACT  (A  brief  and  factual  summary  of  the  document.  It  may  also  appear  elsewhere  in  the  body  of  the  document  itself.  It  is  highly 
desirable  that  the  abstract  of  classified  documents  be  unclassified.  Each  paragraph  of  the  abstract  shall  begin  with  an  indication  of  the 
security  classification  of  the  information  in  the  paragraph  (unless  the  document  itself  is  unclassified)  represented  as  (S),  (C),  or  (U).  It  is 
not  necessary  to  include  here  abstracts  in  both  official  languages  unless  the  text  is  bilingual.) 

The  micro-Doppler  (m-D)  effect  appears  in  the  synthetic  aperture  radar  (SAR)/inverse  SAR 
(ISAR)  image  of  a  target  whenever  the  target  has  one  or  more  rotating  or  vibrating  parts.  M-D 
effect  introduces  distortion  in  the  SAR/ISAR  images.  On  the  other  hand,  m-D  effect  also  carries 
information  about  the  features  of  moving  parts  of  a  stationary  or  moving  target  that  can  be  used 
for  target  identification  purpose.  Based  on  L-statistics,  spectrogram,  and  inverse  Radon  trans¬ 
form,  a  radar  data  exploitation  Matlab  toolbox  was  develped  for  classifying  air,  land  and  ocean 
targets,  which  helps  Geospatial  Intelligence  (GEOINT)  support  to  Canadian  Forces.  This  toolbox 
will  scan  and  search  any  M-D  activity  in  the  desired  area  of  the  SAR/ISAR  image.  The  toolbox 
has  three  options;  1)  Focus  the  stationary  target  if  there  is  no  m-D,  2)  Focus  the  target  after 
removing  the  m-D,  and  3)  extract  the  m-D  and  determine  the  motion  parameters.  The  toolbox 
was  tested  and  validated  against  a  variety  of  simulated  and  measured  targets. 


14.  KEYWORDS,  DESCRIPTORS  or  IDENTIFIERS  (Technically  meaningful  terms  or  short  phrases  that  characterize  a  document  and  could 
be  helpful  in  cataloguing  the  document.  They  should  be  selected  so  that  no  security  classification  is  required.  Identifiers,  such  as 
equipment  model  designation,  trade  name,  military  project  code  name,  geographic  location  may  also  be  included.  If  possible  keywords 
should  be  selected  from  a  published  thesaurus,  e.g.  Thesaurus  of  Engineering  and  Scientific  Terms  (TEST)  and  that  thesaurus  identified. 
If  it  is  not  possible  to  select  indexing  terms  which  are  Unclassified,  the  classification  of  each  should  be  indicated  as  with  the  title.) 

Micro-Doppler;  Fourier  Transform;  Time-Frequency  Analysis;  Short-Time  Fourier  Transform; 
Radon  Transform;  L-Statistics;  Geospatial  Intelligence;  Radar  Detection;  Classification;  SAR; 
ISAR 
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